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ABSTRACT 

We study the distribution of cold dark matter (CDM) in cosmological simulations 
from the FIRE (Feedback In Realistic Environments) project, for M* ~ 10 4-11 M 0 
galaxies in Mh ~ 10 9-12 M 0 halos. FIRE incorporates explicit stellar feedback in the 
multi-phase ISM, with energetics from stellar population models. We find that stellar 
feedback, without “fine-tuned” parameters, greatly alleviates small-scale problems in 
CDM. Feedback causes bursts of star formation and outflows, altering the DM distri¬ 
bution. As a result, the inner slope of the DM halo profile (a) shows a strong mass 
dependence: profiles are shallow at Mh ~ 10 10 — 10 11 M 0 and steepen at higher/lower 
masses. The resulting core sizes and slopes are consistent with observations. This is 
broadly consistent with previous work using simpler feedback schemes, but we find 
steeper mass dependence of a, and relatively late growth of cores. Because the star 
formation efficiency M*/Mh is strongly halo mass dependent, a rapid change in a 
occurs around Mh ~ 10 10 M 0 (M* ~ 10 6 — 1O 7 M 0 ), as sufficient feedback energy 
becomes available to perturb the DM. Large cores are not established during the pe¬ 
riod of rapid growth of halos because of ongoing DM mass accumulation. Instead, 
cores require several bursts of star formation after the rapid buildup has completed. 
Stellar feedback dramatically reduces circular velocities in the inner kpc of massive 
dwarfs; this could be sufficient to explain the “Too Big To Fail” problem without invok¬ 
ing non-standard DM. Finally, feedback and baryonic contraction in Milky Way-mass 
halos produce DM profiles slightly shallower than the Navarro-Frenk-White profile, 
consistent with the normalization of the observed Tully-Fisher relation. 

Key words: galaxies:evolution — galaxies:halos — galaxies:kinematics and dynamics 
— galaxies:structure — dark matter 


1 INTRODUCTION 


Cold Dark Matter with a cosmological constant (ACDM) is 
a successful cosmological model that can simultaneously ex¬ 
plain large scale fluctuations in the cosmic microwave back¬ 
ground and the large-scale structure of the universe that 
forms out of these fluctuations at much later time ( |Spergel| 
et al.|2007; Springel et al.|2005l. However, on much smaller 
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scales, within dark matter halos that host observed galaxies, 
there are indications that the distribution of dark matter is 
inconsistent with the simplest prediction of the cold dark 
matter paradigm. The most obvious and most studied dis¬ 
agreement is in density profiles of dark matter halos inferred 
from observations of dwarf and low surface-brightness galax¬ 
ies. While observed slopes are relatively flat (central density 


slope 

a ~ 0, where central density oc r a ) (jSalucci & Burk- 

ert||2000; Swaters et al. 2003; Gentile et al. 2004, Spekkens 

et al. 2005[ Walter et al. 2008; de Blok et al. 2008; Oh et al. 

2011 

) simulated cold dark matter halos are cuspy (a ~ — 1) 
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(Flores fe Primack|1994 Moore|[l994| Navarro et al.|[l997|. matter core. An alternative mechanism was proposed by 


This problem is known as the cusp/core problem. 

To address this problem, various modifications of dark 
matter properties have been proposed to erase the steep cen¬ 
tral regions and produce a core-like density profile. Exam¬ 
ples include warm dark matter (WDM), whose free stream- 


Dunstan et al. 2011 

Lovell et al. 20121 and self-interacting 

dark matter (SIDM) whose interaction can substantially af- 

feet the central density profile of halos (e.g. Yoshida et al. 

2000 

Burkert))2000 

Kochanek & White 2000; Dave et al. 

2001 

Elbert et al. 

2015 

1. There is still no consensus on 


whether these modifications can solve the problem and sat¬ 
isfy all observational constraints: these modifications can 
produce serious problems on their own. For example |Maccio| 
|et~ al. |201 2a[) found that in order to produce dark mat¬ 
ter cores as large as those seen in observed dwarf galax¬ 
ies, the warm dark matter would also prevent the forma¬ 


tion of dwarf galaxies. Simple SIDM models ( Carlson et al. 
1992 Machacek et al.|1993 de Laix et al.|1995 l were shown 

to violate observations of central regions of galaxy clusters 
( Miralda-Escude1|2002 Yoshida et al. 20001 that are found 
to be denser and more elliptical than SIDM would predict. 
However, recent SIDM models that take into account more 
accurate observational constraints and the effects of baryons 
offer promising explanation of the problem without violat¬ 
ing any known observational constraints (Rocha et al.|2013 


Peter et al.|2013 Kaplinghat et al.|2014 Elbert et al.|2 015). 
Initial problems with the SIDM have motivated more 


complex models such as velocity-dependent SIDM (Yoshida 
|et al. | [2000[ |Loeb fe Weiner]|2011| |Maccid et al. 2012a|. 


However, SIDM with a simple power-law velocity depen¬ 
dence will not be able to create a core and, at the same 
time, produce stable halos of dwarf galaxies over a Hubble 
time ( Gnedin fe Qstriker|2001 1. Loeb & Weiner (20111 pro¬ 
posed SIDM with a Yukawa potential, which has a nontrivial 
velocity-dependence that is effective at producing cores in 
dwarf galaxies without adverse effects on clusters of galax¬ 
ies. Cosmological simulation of a Milky Way-mass halo with 
this SIDM showed that realistic cores can be formed in sub¬ 


halos expected to host dwarf galaxies (Vogelsberger et al. 


2012a). However, this model requires more free parameters 


for the velocity dependence and it is not yet known whether 
it can reproduce the correct halo abundance and mass dis¬ 
tribution. 

Before concluding that simple cold dark matter mod¬ 
els must be modified,, we must also examine the effects of 
baryons on the distribution of dark matter within halos. 
Baryons are not only what is actually observed in galaxies, 
but baryonic effects at the halo center can, in principle, also 


play a role in shaping the dark matter profiles (Blumenthal 

et al 

1986| |Navarro et al. 1996; El-Zant et al. 2001; Gnedin 

et al 

2004; ^eaTnr Gilmore 2005; Governato et al. 2010 

Peiiarrubia et al. 2012; Governato et al. (2012, Pontzen & 

Governato 2012; Maccio et al. 2012b; Teyssier et al. 2013; 

Di Cintio et al.|2014 

Pontzen & Governato|2014). 


'Javarro et al. 

(1996) and Read & Gilmore (2005) 


used N-body simulations to model a sudden removal of 
a large baryonic component via supernova-driven winds 
(represented as a change in the external potential) in a 
dwarf galaxy with an initially cuspy dark matter halo. They 
showed that such mass removal leads to formation of a dark 


Mashchenko et al. (20061, who showed that bulk motion of 


gas within forming galaxies leads to significant gravitational 
potential changes which can also redistribute dark matter 
and reduce its central density. Dynamical effects between 
baryons and dark matter during the halo formation were 
also suggested as a mechanism that could modify dark mat¬ 


ter density profiles (e.g. El-Zant et al. 2001 Tonini et al. 
2006 Romano-Dfaz et al.|2008 Del Popolo|2009 |. 

Recently Governato et al. f( |2010| |2012[ ) used cosmologi¬ 
cal “zoom-in simulations” with baryons, cooling, star forma¬ 
tion and supernovae feedback to show that outflow episodes 
in dwarf galaxies can turn the central dark matter cusps into 
cores. Strong supernova-driven outflows from clustered star 
formation in the inhomogeneous interstellar medium (ISM) 
resulted in a decrease of the dark matter density within cen¬ 
tral kiloparsec to less than half of what it would otherwise 
be in a halo of this mass (md m ~ 10 10 Mq). Pontzen & Gov 


ernato 


(2012) further clarified this density flattening mech¬ 


anism: a quick change in gravitational potential due to gas 
outflow can effectively inject energy into dark matter or¬ 
bits and (typically after many outflow episodes) flatten the 
central dark matter profile. They showed that the repeated 
changes of gravitational potential on timescales shorter than 
tdyn during 2 < z < A can significantly flatten cuspy dark 
matter profiles. 


Brook et al. (2012) showed that a large fraction of the 


gas that is expelled returns via a large-scale galactic fountain 


(see also Oppenheimer et al. 20101 to form stars at later 
times: this greatly increases the chance of outflows from the 
inner regions and further helps in core formation. Other, 
non-cosmological simulations with strong SNe feedback also 
showed that it is possible to form cores in dwarf galaxies 
owing to bursty star formation that removes large quantities 
of gas during bursts (Teyssier et al. [2013 1. 

On the other hand, idealized simulations by|Gnedin fc 


Zhao 


et al. 


(2002 


) Ogiya fe Mori| (|2011[), and | Garrison-Kimm el 


120131 focused on time evolution, supernovae energy 


requirements and mass ejection frequency in idealized mod¬ 
els and argued that SNe driven feedback is not efficient 
enough to form cores at the observed level. However, re¬ 
cent cosmological simulations have had more success. |Madau| 
et al.|p014|) suggest that earlier mass removal can lower the 
energy requirements for core formation and |Onorbe et al.| 
(20151 showed that late star formation, after the early epoch 


of cusp building, is particularly efficient at utilizing stellar 
feedback to remove dark matter. 

There are two other related “problems” with structural 
properties of CDM halos. One is the lack of very steep 
central profiles in relatively massive disk galaxies. Cuspy 


NFW profiles (Navarro et al. 19971 are expected to be even 


steeper within baryon dominated galaxies owing to the con¬ 
traction of dark matter caused by the central concentra¬ 


tion of baryons (Blumenthal et al. 1986). The distribution 


of matter in galaxies effectively determines the Tully-Fisher 
relation (e.g. Tully fe Fisher|1976 Dutton et al.|2007 1. How¬ 
ever, given the observed distribution of baryons, contraction 
of an NFW halo would result in circular velocities too high 
at a given luminosity/mass. This motivated several authors 
to suggest that dark matter does not undergo contraction 
or is perhaps even expanded from the original cuspy NFW 
profile (Dutton et al. 2007). While this could be interpreted 
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as suggesting a problem with the currently favored CDM 
model, stellar feedback is also able to effectively “expand” 
the dark matter distribution even in Milky-Way mass halos 
( Maccid et al.|2012b l. 

The second problem is the so called “Too Big To Fail" 


problem (Boylan-Kolchin et al. 2011). In the Milky Way, 


satellites have significantly lower dark matter densities in 
the inner few hundred parsecs than the corresponding sub¬ 
halos in CDM only simulations without baryons. Alterna¬ 
tively, massive sub-halos whose inner densities are high, 
never formed galaxies. Similar problems also exist in other 


dwarf galaxies in the Local Group (Garrison-Kimmel et al. 
|2014| [Papaster gis et al.|2015| ). There are hints that feedback 
can help solve this issue along with the cusp/core problem 


I Madau et al.|[2014 Brooks fe Zolotov||2014| Onorbe et al. 
[201^1, although proper statistics are still lacking. 

It is becoming clear that the bursty nature of stellar 
feedback in galaxies can modify the inner regions of dark 
matter halos. However, in general, most simulations used 
to study this problem so far used crude and often unphys¬ 
ical implementations of stellar feedback. One might worry 
that this could impact the effect of stellar feedback. Most of 
the cosmological simulations used to address the cusp/core 
issue simply turn off cooling from the gas heated by su- 


pernovae ejecta until such gas escapes galaxies ( 

Governato 

et al. 2010. 2012 J |Pontzen & Governato 2012;|Maccio et al. 

2012b 

Di Cintio et al. 

2014). The delayed cooling is un- 


physically long and results from a misinterpretation of the 
standard supernova remnant results (Martizzi et al.| |2015[ ). 
In addition, most simulations include only supernovae feed¬ 
back while other stellar feedback mechanisms are ignored or 
implemented crudely: e.g. radiation pressure, cosmic rays, 
and photo-heating are often approximated with pure ther¬ 
mal energy input and additional freely-adjustable parame¬ 
ters ( Maccid et al.|2012b Di Cintio et al.|2014 l. 

Furthermore, the particle mass resolution used in some 
previous studies was insufficient to properly resolve the ob¬ 
served core sizes. Low resolution may hinder the investiga¬ 
tion of central density profiles on small scales owing to two- 


body relaxation effects (Power et al. 2003). In Appendix 


[A] we show the relation between the convergence radius and 
particle mass that should be used to estimate resolved scales 
in different simulations. 

To study the cusp/core problem in a complex cosmologi¬ 
cal environment we use simulations fro m the Feedba ck in Re¬ 
alistic Environments (FIRE) project (Hopkins et al. .2014 
(H14), supplemented by four new dwarf galaxy simulations. 
Our simulations use physically motivated stellar feedback, in 
which energy and momentum input are based on stellar pop¬ 
ulation synthesis models alone, with no adjusted parameters. 
In additions to supernovae energy and momentum we in¬ 
clude radiation pressure, stellar winds, photo-ionization and 
photo-electric heating processes. In H14 we show that the 
M*-Mh relation in FIRE is in reasonable agreement with ob¬ 
servations, for galaxies residing in halo masses M < 10 12 Mq. 
This result is sensitive to the feedback physics: simulations 
with supernovae alone fail to reproduce the correct rela¬ 
tion, unless additional feedback processes are also incorpo¬ 
rated. Overall, stellar feedback in FIRE simulations results 


Project website: http://fire.northwestern.edu/ 


in bursty star formation histories followed by strong out¬ 
flow episodes ( Muratov et al.| [2015l that can affect matter 
distribution within galaxies. 

Our simulations are among the highest resolution cos¬ 
mological “zoom-in" simulations to date evolved down to 
z = 0 with full baryonic physics. In addition to the ad¬ 
vantages in implemented physics and resolution, we adopt 


the P-SPH “pressure-entropy” formulation of SPH (Hopkins 


20131, which includes a large number of numerical improve¬ 


ments relative to previous SPH studies, which together sig¬ 
nificantly improve the treatment of cooling and multi-phase 
fluid mixing, and reduce the well-known discrepancies be¬ 
tween SPH and grid-based codes. 

In this paper we study halos with masses 10 u < 

Mj,/Mq < 10 12 with full feedback and their dark matter 
only analogs, which enables us to directly compare their 
dark matter distributions. We find results in broad agree¬ 
ment with previous work, but with some important differ¬ 
ences. We find that stellar feedback affects all of the sys¬ 
tems we study but large cores develop only in the halo mass 
range of ~ 10 10 Mq to a few x 10 11 Mg. Furthermore we show 
that cores change over time, and that progenitors of massive 
galaxies once had more prominent cores. We demonstrate 
how bursty star formation and related feedback correlate 
with changes in dark matter halos and show that feedback ef¬ 
fectively cancels the effects of adiabatic contraction. Finally 
we discuss consequences of our results for the cusp/core is¬ 
sue, the Tully-Fisher relation, the “Too Big To Fail" problem 
and indirect dark matter detection. 

We find encouraging trends that have the potential to 
solve most of the apparent small scale problems of the CDM 
paradigm. 

The paper is organized as follows: includes a brief 

description of the code and implemented stellar feedback as 
well as the set-up of the simulations. In © we show the 
dark matter density profiles and their time evolution. In 
^4] we study the effects of stellar feedback on the expected 
contraction of dark matter and on the Tully-Fisher relation. 
In ij5] we compare our results with previous work, discuss the 
implications and propose directions for further investigation. 


2 SIMULATIONS 
2.1 Simulation Code 

The simulations in this work were run with the newly devel¬ 


oped GIZMO code (Hopkins 2014) in a fully conservative, 


pressure-entropy based smoothed particle hydrodynamic (P- 
SPH) mode (H opkins|2013 1. P-SPH eliminates artificial sur¬ 
face tension at contact discontinuities that affects traditional 
density based SP H (|Agertz et al.||2007| |Saitoh fe Makino| 
2013 Sijacki et al.||20f2 l. We use the artificial viscosity al¬ 
gorithm with a switch from Cullen & Dehnen (20101 which 


reduces viscosity to close to zero away from shocks and en¬ 
ables accurate shock capturing. The same higher-order dis¬ 
sipation switch is used to trigger entropy mixing at the ker¬ 
nel scale following P rice| (|2008[). Time-ste pping is controlled 
by the limiter from Durier & Dalla Vecchiaj (20121, which 


limits the difference in time-steps between neighboring par¬ 
ticles, further reducing numerical errors. The gravity solver 
of the GIZMO code is an updated version of the PM+Tree 
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4 T. K. Chan et al. 


algorithm from Gadget-3 (Springcl 20051 and uses fully con¬ 
servative adaptive gravitational softening for gas ( Price fe 
Monaghan 2007). GIZMO’s softening kernel represents the 
exact solution of the particle mass distributed over the SPH 
smoothing kernel (B arnes|2012 1. 

The code performs well on standard strong shock, 
Kelvin-Helmholtz and Rayleigh-Taylor instabilities, and 


subsonic turbulence tests (for more details see Hopkins 


20131. In cosmological “zoom-in" simulations of a Milky Way 
size halo without outflows, the code eliminates most of the 
artificial fragmentation of halo gas seen in traditional SPH 


simulations (Kaufmann et al. 2006 Sommer-Larsen 2006 


Keres & Hernquist 2009) and increases cooling from the hot 


halo gas at late times (e.g. jKeres et al. 20051, when compared 
to the classical SPH (Keres et al., in preparation). Overall, 
the resulting halo gas properties are similar to the results 


from adaptive-mesh and moving-mesh simulations (Agertz 


et al.|2009| [Keres et al.|2012[|Vogelsberger et al.|2012b[|N eI- 

son et al.||2013 1. 


2.2 Baryonic Physics 

Our simulations incorporate cooling, star formation and 
physical stellar feedback processes that are observed to be 
relevant in the inter-stellar medium. Here we briefly review 
these components, for detailed description please see H14. 

Gas follows an ionized+atomic+molecular cooling curve 
from 10 — 10 10 K, including metallicity-dependent fine- 
structure and molecular cooling at low temperatures, and 
at high-temperatures (> 10 4 K) metal-line cooling followed 
species-by-species for 11 separately tracked species. At all 
times, we tabulate the appropriate ionization states and 
cooling rates from a compilation of CLOUDY runs, includ¬ 
ing the effect of the photo-ionizing background. We use 
global ultraviolet background model from [Fauchcr-Giguere| 


et al. (2009) that heats and ionizes the gas in an ionization 


equilibrium approximation. We apply on-the-fly ionization 
corrections in denser gas to account for the self-shielding 
based on the local Jeans-length approximation (to deter¬ 
mine the surface density), which provides an excellent match 


to a full ionization radiative transfer calculation (Faucher- 


Giguere et al.|[2010| |Rahmati et al.||2013| |Faucher-Giguere 
et al.|20~ I. 


Star formation is allowed only in dense, molecular, self- 
gravitating regions above n > 10—100 cm -3 . This threshold 
is much higher than that adopted in most “zoom-in" simu¬ 
lations of galaxy formation (the high value allows us to cap¬ 


ture highly clustered star formation). We follow Krumholz 
|fe Gnedln| ( |2011| ) to calculate the molecular fraction /h 2 in 
dense gas as a function of local column density and metallic- 
ity. We allow SF only from the locally self-gravitating molec¬ 
ular gas using the efficiency of 100% per free fall time (the 
actual SF efficiency is feedback regulated). 

Our stellar feedback model includes a comprehensive 
set of physical mechanisms: radiation pressure, supernovae 
(with appropriate momentum and thermal energy input), 
stellar winds, photo-ionization and photo-electric heating 
as described in H14. We do not tune any feedback model 
parameters but instead directly use the energy, momen¬ 
tum, mass and metal return based on the output of the 
STARBURST99 stellar population synthesis model ( |Lei-| 
therer et al. 19991. Our feedback model is implemented 


within the densest interstellar-medium material, yet we do 
not resort to turning off cooling of supernova heated gas at 
any time. 

2.3 Initial conditions and zoom-in method 

We adopted a ‘standard’ flat ACDM cosmology with flo ~ 
0.27, A « 0.73, flj, « 0.045, h ~ 0.7 for all runs. In order to 
reach the high-resolution necessary to resolve a multi-phase 
ISM and to properly incorporate our feedback model we 
use the “zoom-in” technique. This places maximum baryonic 
and dark matter resolution around the halo of interest in 
a lower resolution, collisionless box (Porter 1985 Katz & 
Whi te|1993 1. 


We consider halos with mass from 10 9 to 10 12 Mq at 
z = 0 from the FIRE project (Hopkins et al. 2014). Ini¬ 
tial conditions of those halos are listed in Tab. [7] The sim¬ 
ulations m09 and mlO are constructed using the meth¬ 


ods from Onorbe et al. (2014); they are isolated dwarfs. 


Simulations mil, ml2q and ml2i are chosen to match a 


subset of initial conditions from the AGORA project (Kim 


et al. 2014I while ml2v uses initial conditions from Faucher 
Giguere & Keres (2011) (higher resolution versions of the 
run first presented in |Keres fc Hernquist 20091. In addi¬ 
tion, we re-simulated all of these initial conditions using 
dark matter only N-body simulations with the same to 
have a matched set of simulations with and without baryonic 
physics for a direct comparison. 

To improve halo mass coverage in the regime where 
cores are prominent, Mh ~ 10 1() — 10 41 Mq, we have also sim¬ 
ulated additional four halos with the same FIRE code, also 
listed in Tabled These additional simulations are first time 
presented in this work. Initial conditions were generated us¬ 
ing the MUSIC code ( Hahn fe Abel|2f)TT I. We randomly se¬ 
lected halos with small Lagrangian regions for resimulation 
from a 40 h -1 Mpc box. All particles within 3R VIT at z = 0 
are enclosed by the Lagrangian region to reduce contami¬ 
nation (Onor be et al.|2 014l. Halos ml0hl297, ml0hll46 
and ml0h573 have masses between mlO and mil, while 
mllh383 is slightly more massive. They are all isolated 
field dwarfs. 


2.4 Convergence Radius 


We adopted the method described in Power et al. (2003) to 


calculate a conservative limit for the convergence radius of 
the dark matter profiles in the N-body only simulations. 
They found that effective resolution is related to radius 
where the two-body relaxation time, t re iax, becomes shorter 
than the age of the universe to- They verified this with N- 
body simulations and found out that for this particular prob¬ 
lem, well resolved regions of halos require t re i ax > 0.6to- At 
smaller radii, even if the dynamics are locally well resolved, 
small N-body effects can, over a long integration time, arti¬ 
ficially turn cusps into cores. Given the enclosed number of 
particles N and the average density of the enclosed region p 
one can show that: 

- 1/2 

(1) 


Uelax(r) _ C200 N f p \ 

to 8 In N V Pcrit ) 


where p cr it is the critical density. We define rp ow as the 
smallest radius that fulfills t re iax > 0.6£o for the dark mat- 
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Name 

M h u 


Rvir 

Rl/2 

m b 

Cb 

m dm 

£dm 


[M e ] 

[M 0 j 

[kpc] 

[kpc] 

\Mq\ 

[pc] 

[M 0 ] 

[pc] 

m09 

2.6e9 

4.6e4 

36 

0.49 

2.6e2 

2.0 

1.3e3 

29 

mlO 

7.9e9 

2.3e6 

50 

0.51 

2.6e2 

2.9 

1.3e3 

29 

mil 

1.4ell 

2.4e9 

1.4e2 

6.9 

7.1e3 

7.1 

3.5e4 

71 

ml2v 

6.3ell 

2.9el0 

2.2e2 

1.8 

3.9e4 

21 

2.0e5 

2.9e2 

ml2i 

l.lel2 

6.1el0 

2.7e2 

4.3 

5.7e4 

20 

2.8e5 

1.4e2 

ml2q 

1.2el2 

2.1el0 

2.8e2 

3.6 

7.1e3 

20 

2.8e5 

2.1e2 

dm09 

3.3e9 

- 

39 

- 

- 

- 

1.6e3 

29 

dmlO 

9.3e9 

- 

54 

- 

- 

- 

1.6e3 

29 

dmll 

1.6ell 

- 

1.4e2 

- 

- 

- 

4.3e4 

71 

dml2v 

7.7ell 

- 

2.4e2 

- 

- 

- 

2.4e5 

2.9e2 

dml2i 

l.lel2 

- 

2.9e2 

- 

- 

- 

3.4e5 

1.4e2 

dml2q 

1.4el2 

- 

3.0e2 

- 

- 

- 

3.4e5 

2.1e2 

ml0hl297 

1.3el0 

1.7e7 

62 

1.8 

1.5e3 

4.3 

7.3e3 

43 

ml0hll46 

1.6el0 

7.9e7 

65 

2.5 

1.5e3 

4.3 

7.3e3 

43 

ml0h573 

4.0el0 

3.2e8 

88 

3.4 

1.5e3 

10 

7.3e3 

1.0e2 

mllh383 

1.6ell 

4.0e9 

1.4e2 

7.2 

1.2e4 

10 

5.9e4 

1.0e2 

dml0hl297 

1.6el0 

- 

66 

- 

- 

- 

8.8e3 

43 

dml0hll46 

1.8el0 

- 

68 

- 

- 

- 

8.8e3 

43 

dml0h573 

4.2el0 

- 

91 

- 

- 

- 

8.8e3 

le2 

dmllh383 

1.7ell 

- 

1.4e2 

- 

- 

- 

7.1e4 

1.0e2 


Table 1. Simulation details. and M* are the total mass and stellar mass of the largest halo in the simulation at z = 0. R ±/2 is the 
radius of the region where half of the stellar mass is enclosed, raj-, is the mass of a gas particle in the simulation; ra^m is the mass of a 
dark matter particle in the simulation. €5 is the minimum gravitational smoothing length of gas; €d m is the Plummer equivalent 
gravitational smoothing length of dark matter. The simulation name convention is as follows, “mXX” refers to the halo mass 
- 10 xx M o . d are the corresponding dark matter only simulations, e.g. dm09 corresponds to m09. DM only and hydro dynamical 
simulations have the same initial conditions, except that the gas particles are absorbed into the dark matter particles in DM only 
simulations. m09, mlO, mil, ml2v, ml2i and ml2q are from H14, whereas ml0hl297, ml0hll46, ml0h573 and mllh383 are 
new simulations first presented in this work. 


ter only simulations and use rp ow to conservatively estimate 
where DM profiles are converged. We show the minimum 
particle mass required for converged profiles at 0.3 —3%7? v ir 
in Appendix [A] and discuss implications and limitations of 
this convergence criterion in § |5.4| 


2.5 Halo Finding 


We identify halos and estimate their masses and radii us¬ 


ing the Amiga Halo Finder (AHF) (Knollmann & Knebe 
20090 AHF uses an adaptive mesh refinement hierarchy 
to locate the prospective halo center (Knebe et al. 20011. 


We use the Bryan & Norman (19981 formulae to determine 
the virial over-density and virial radius, R v i r , and quote the 
halo mass as the mass enclosed within the virial radius. We 
follow the main progenitor of a halo using the merger tree 
code included in AHF and check the growth history of each 
individual halo, making sure that we follow the same main 
progenitor. We use the histories of the main progenitors to 
study the time evolution of the density profile. Occasionally, 
AHF misidentifies the main progenitor of the halo, so some¬ 
times 7? v ir temporarily decreases. We avoid this problem by 
searching for another halo with larger R v i r and its center 
within 50 kpc of the center of the main progenitor in the 
previous snapshot. 

During ongoing mergers and for galaxies with large 
clumps of stars and gas, AHF adaptive centering on the 


2 http://popia.ft.uam.es/AHF/Download.html 


highest overall density might quickly change over time and 
might not center on the stars. This is especially important in 
halos with shallow dark matter profiles and relatively shal¬ 
low distributions of stars. To avoid this issue and have a con¬ 
sistent centering on the stellar component, we use two-step 
procedure to identify the center. First we use AHF to de¬ 
fine Rvir and an approximate center. Within 0.U? v ir around 
this center we place the stellar mass volume density on a 
grid. We search for progressively higher overdensities which 
we enclose in an iso-density ellipsoid. Once this ellipsoid in¬ 
cludes less than one quarter of the total galaxy stellar mass 
we stop the procedure. The center of the ellipsoid is our new 
halo center. This ensures consistent centering of all profiles 
on the galactic stellar distribution. We tested this proce¬ 
dure and found that our newly defined center is closer to 
the DM density peak than the original AHF for cases with 
shallow central cores and shows less stochastic variations of 
the central density slope over time. 


3 DARK MATTER PROFILES AND THEIR 
TIME EVOLUTION 

The focus of this paper is the effect of stellar feedback on the 
dark matter distribution in our simulated halos. There are 
three major areas we explore: the relation between the halo 
mass and the central dark matter profile, the time evolution 
of the inner dark matter density, and the changes to galaxy 
structure caused by the effects of stellar feedback. 


© 0000 RAS, MNRAS 000, 000-000 



















6 T. K. Chan et al. 








Figure 1 . Dark matter density profiles of halos at z = 0. Black dashed lines represent collisionless dark matter only simulations; red 
solid lines represent simulations with baryons and stellar feedback. The Power radius rp ow , within which N-body relaxation effects can 
become important, is shown with vertical black dashed lines. The halo masses are shown in the brackets. Baryonic feedback reduces the 
central DM density, especially at around M\ t ~ 10 11 Mg. 


3.1 Dark matter profiles 

Figure[l]shows spherically averaged dark matter density pro¬ 
files of six simulated halos at z = 0. We focus on the inner re¬ 
gions of halos, 0.002-0.2f? v ir where galaxies reside and where 
effects of feedback are expected to be measurable. We show 
both the profiles from DM only simulations as well as DM 
profiles from simulations with baryonic physics. DM only 
profiles are re-normalized to account for the lower global 
Odm in simulations with baryons. Effects of the baryonic 
physics are visible in most halos to a different degree and 
are the largest in mil and ml2v simulations. In m09 the 
DM density profile is almost the same in simulations with 
and without baryons, while in mlO, a small, resolved core 
forms in the central region. The density profile in mil has 
the largest core and the lowest central DM density. In ml2v, 


central density starts increasing again and the relative core 
size decreases, but differences in profiles are present all the 
way to several percent of R v i r . In the two most massive halos 
we analyze, ml2i and ml2q, differences in central region 
are even smaller, although profiles are still shallower than 
what is expected based on N-body simulations. However in 
§ |3.5[ we show that the effect of feedback in these halos is 
significant and largely cancels out the gravitational influ¬ 
ence of baryons that is expected to steepen the profile seen 
in dark matter only case. We note that all of the plotted 
range is resolved with many gravitational softenings of the 
dark matter particles. We also show the more conservative 
Power convergence criterion which is typically a fraction of 
a percent of R v n for all of our halos. 

We quantify the effect of feedback on dark matter distri¬ 
bution using two parameters that are frequently estimated 
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Figure 2. Dark matter density profiles of halos at z = 0 from hy¬ 
drodynamic simulations (blue solid circles), fits with the pseudo- 
isothermal sphere (black solid line) and fit with a power-law model 
(oc r“) at 1-2 %R vlr (green dashed line). The black vertical line 
shows the convergent radius according to the Power criterion. The 
pseudo-isothermal sphere is a good fit to the central regions of the 
simulated halos and provides a good estimate of the core sizes. 


from observations: the central slope a of the dark matter 
density profile (pDM oc r a ) and the core radius r core of the 
pseudo-isothermal sphere (see Eq. [2]). Examples of the fits 
are shown in Figure [2] 


3.2 Inner slopes of dark matter halo profiles 


We estimate the slope a of the dark matter density profile 
by fitting a power law relation pDM oc r a in the 1 — 2%R V1I 
interval. This range is well resolved for all of our main halos 
at z=0 and it is physically meaningful as it shows the rel¬ 
ative profile change at a fixed fraction of the halo size. For 
dwarf galaxies, this is close to the region where observations 
indicate shallow and core-like profiles in low-mass galaxies 
(typically measured at a few hundred pc; see Oh et al.|2011 


|2015| |Walter et al.|2 0081 Hu nter et al.||2012 1. We also show 

slopes at 0.5— l%7? v i r for comparison. Example of the fit for 
a are shown in Figure [2] We have also measured a in the 
fitting ranges of 0.3 — 0.7kpc and 1 — 2kpc. In Appendix [B] 
we discuss limitations of these alternative fitting ranges and 
show that general trends of a with halo mass are similar to 
our default fitting choice. 

Figure [3] shows a as a function of the halo virial mass, 
Mh, at different redshifts. We only show main halos with 
more than 10 5 DM particles and remove all sub-halos and 
halos with more than 1% contamination by more massive 
DM particles within the inner 0.1-R v ir[^] The hollow circles 
represent a whose fitting range contain regions smaller than 
0.5 rp ow and/or larger than 1/3 r s , where r s is the scale 
radius of the NFW profiles. Q 

We focus on z ^ 2 when profiles of halos start to stabi¬ 
lize as rapid halo growth subsides. At z = 0, the simulated 
halos show a clear tendency to form shallow central profiles 
at Mh ~ 10 10 — 10 11 Mq. All of the profiles in this range 
are significantly shallower than the NFW profile. More ac¬ 
curate estimate of the halo mass and stellar mass ranges 
where feedback flattens central slopes will require a larger 
number of simulations as our statistic are currently limited. 
When profiles are measured at even smaller radii, 0.5-1% 
of AMr, profiles are typically even more shallow. At z=2 we 
see that the scaling with mass shows much larger dispersion, 
which owes to very bursty star formation and central halo 
regions that are just coming out of the fast growth stage. 
We later show that in intermediate mass halos at a fixed 
physical radius, DM profiles get shallower with time. 

It is interesting to notice that low mass dwarfs with 
Mh -C 10 10 Mq do not develop density cores even at 1% of 
R v ir (which is typically only several hundreds of parsecs). As 
we discuss later, only a small fraction of baryons are con¬ 
verted to stars in these halos, owing to efficient feedback and 
effects of the UV background. The energy available from a 
small number of SNe is not sufficient to dramatically mod¬ 
ify the dark matter distribution. Around Mh = 10 10 Mq, 


3 Note that in most runs we use a “padding region” around our 
zoom in region where mass resolution is lower only by factor of 
8 . Mild contamination with such particles can sometimes occur 
within R v ; r but typically has no consequence on the evolution of 
halo gas or central slopes. 

4 We estimate r s from the concentration c = R v i T /r s at a given 
mass, from the concentration-mass relation of|Dutton & Maccio| 
12014 i. 
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Figure 3. Slopes of dark matter density a as a function of halo 
mass at different redshifts. The green dotted and red dashed lines 
represent the expected a for the NFW profiles measured at 1- 
2% and 0.5-l%_R v ir respectively. The concentration of the NFW 
profiles is evolved with redshift as in [Dutton fe Maccio| (2014). 
Filled circles represent a for simulated halos in which the fitting 
range is larger than 0.5 rp ow and smaller than one third of r s . 
Hollow circles represent the slopes in halos in which at least one of 
these criteria is not satisfied (see the main text for other selection 
criteria) At ~ 10 10 — 1O 11 M0, baryonic effects lead to profiles 
significantly shallower than the corresponding NFW profiles from 
N-body simulations. 
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Figure 4. ( Upper) Relation between a and ratio between M* and 
Mi,. The blue dashed line show the fit from Di Cintio et al. j ([2014) 
whereas the yellow dotted line shows the same fitting for our data. 
(Lower) Relation between a and M*. The symbols are explained 
in Figure [3] The DM profiles near halo centers are cuspy at the 
lowest and highest masses, and shallowest at M* ~ 10® — 10 9 Mq 
and ~ 0.01. 


the slope of the inner density profile increases rapidly with 
mass, indicating the development of DM cores. This seems 
to be a “threshold” halo mass needed to develop large cores. 
As discussed in Onorbe et al. (20151, small differences in 
star formation histories in halos close to this threshold can 
result in the substantial difference in central slopes of the 
dark matter distribution. 

Finally, in halos with mass comparable to the Milky 
Way (ml2v and ml2i) profiles steepen again and are only 
slightly shallower than NFW. These halos have deep poten¬ 
tial wells that can retain a large fraction of available baryons 
and convert them into stars. Baryons are actually expected 
to steepen the DM profiles to a < —1 owing to adiabatic 
contraction of dark matter. However, bursty feedback largely 
cancels and in some cases even overcomes this expected ef¬ 
fect of contraction, resulting in slopes a > — 1. The inter¬ 
play between baryonic contraction and stellar feedback will 
be discussed in Section f3.51 

Figure [4] shows the scaling of profile slope with galaxy 
stellar mass (lower panel) and M*/Mh (upper panel). In 
terms of stellar mass, feedback significantly modifies DM 
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slopes in the M* ~ 10 7 — 5 x 10 9 Mq range, with a fast tran¬ 
sition from cusps to cores occurring at a few x 10 6 — 10 7 Mq. 
Overall trends of a with Mh in Fig ure [3] are similar to the 
result of Governato et al. (20121 and Di Cintio et al. (20141 


(see also a recently submitted work by Toilet et al. 


20151. 


However, we stress that both of these simulations simply 
suppress cooling in dense gas after supernovae explosions 
rather than explicitly treat most of the feedback processes 
around young stars. Furthermore the spatial and mass res¬ 
olution is typically better in our simulations, by about a 
factor of ten in mass. This leads to some differences in the 
slopes of dark matter halos that are illustrated in the upper 
panel of Figure [4] In general, profile slopes increase faster 
with M* / M h compared to the previous “subgrid” models, 
suggesting faster transition from cusps to cores. We cau¬ 
tion that a small number of simulated halos in both samples 
could also be responsible for some of the differences. In addi¬ 
tion we find that the central slope relation is different for the 
inner 1% or inner 2% of _R v ir, which means that the fitting 
formula in Di Cintio et al. (20141 is not generally applicable. 
We compare our result with observations in §4.3| 


3.3 Core radii 


In addition to the inner slope, we also examine another pa¬ 
rameter, the core radius r co re of the halo. We quantify the 
core size using the the pseudo-isothermal sphere fit that is 
frequently used to describe dark matter density profiles (e.g. 


Begeman|1987[|Broeils|1992l|de Blok fe McGaugh|1997HVer-| 

heijen||1997| anoth er popular fit is the Burkert profile, e.g. 

Salucci & Bur kert| |2000). Density profile is given by 


p{r) = Po 


1 + 


( 2 ) 


where po is the central dark matter density. We use Eq. 
[2] to fit the spherically averaged dark matter density pro¬ 
files of the simulated halos. The two free parameters, r CO re 
and po, were determined through a y 2 minimization fit¬ 
ting procedure starting at r = O.lkpc and ending at r = 
min[R v i r , lOOkpc]. Table [ 5 ] lists r core for all halos analyzed 
in this work. Examples of our fits are shown in Figure [2] In 
general fits agree well (to better than few tens of percent 
within 0.1-Rvir) with the DM density profiles for all halos 
with Mh < 10 12 Mq. For ml2q and ml2i pseudo-isothermal 
profiles deviate significantly from the DM distribution. Evi¬ 
dence of cores is present in all halos. The core size is smallest 
(relative to R v i r ) in the m09 run (< 0.5%_R v ir) and largest 
in mllh383 and ml0h573 (> 4%7? v i r ) , where we also find 
the shallowest central slope. Cores of the size of > 0.005i? v ir 
are present even in Milky Way mass halos, albeit with higher 
central DM density and less shallow central slopes than at 
Mh ~ 10 11 Mq. We compare the core radii from our simula¬ 
tions to observations in § |4.3| 


3.4 Time Evolution of a 

Next, we investigate the time evolution of the central slope 
of the dark matter distribution, for five representative ha¬ 
los from our sample. Left panels of Figure [5] show the time 
evolution of a measured at 1 — 2% of the halo virial radius 
at 2 = 0 (Rvo)- For each halo, this radius is kept fixed in 


2 = 0 

r core (kpc) 

T core / Rvir 

m09 

0.17 

0.0048 

mlO 

0.38 

0.0073 

mil 

4.7 

0.034 

ml2v 

1.4 

0.0061 

ml2i 

2.0 

0.0073 

ml2q 

1.2 

0.0043 

ml0hl297 

2.0 

0.032 

ml0hll46 

2.1 

0.033 

ml0h573 

3.6 

0.041 

mllh383 

5.7 

0.041 


Table 2. The core radii of the best fitted pseudo-isothermal 
spheres (Eq. [2j of the simulated halos at z = 0. Large cores of 
3-4% -R v i r form at Mh ~ 10 10 — 10 11 Mq. 


physical units at all times. [^In mlO, a < —1 at all times, 
which is consistent with its relatively small core size, below 
1%-Rvir- In ml0hl297, a is steadily rising from — 1 at z> 1 
to ~ —0.5 at 2 = 0. In mil, a ~ —1 early on, increasing to 
a ~ 0 around z ~ 1 and stays quite flat until late timesQ 
In ml2i and ml2v central DM profiles are flattened to 
a ~ —0.7 at z ~ 1 but steepen afterwards such that central 
dark matter slope is a ~ —1 at the present time. 

The right panels show the time evolution of the en¬ 
closed mass within 0.027? v o and the star formation rate 
within 0.1-RvO averaged over 0.2 Gyr. A dense central 
concentration of dark matter builds up early in mlO, with 
some fluctuations during the bursty star formation epoch at 
2 < z < 4, but as star formation subsides the amount of dark 
matter in the central region remains almost unchanged un¬ 
til the present time. The correlation between star formation 


and strong outflows of gas that follow the burst (Muratov 
et al. 20151 and the decrease of dark matter in the central 


region is clearly visible in the ml0hl297 and mil panels. 
Removal of DM mass occurs after strong bursts of star for¬ 
mation. We examine this more closely in §5.2[ While some 
of the dark matter gets re-accreted, the central concentra¬ 
tion of dark matter remains lower after the burst for at least 
several Gyrs. 


5 We have also analyzed the results within a fixed fraction of 
the time dependent virial radius (R v ir instead of R v o) but the 
correlation between stellar feedback and the enclosed DM mass 
is difficult to interpret because enclosed DM mass increases with 
Rvir • 

6 At late times z < 0.5, this halo undergoes several episodes of 
very fast central slope variations. We examined this system closely 
and found that these are caused by close passages of a substruc¬ 
ture in an ongoing merger. Vertical lines in Figure [5| indicate the 
times of closest passages: they correlate well with temporary drop 
and strong oscillations in central slope. The close passages can af¬ 
fect the accuracy of locating center of the galaxy in AHF, which 
further motivates our two-step center finding procedure described 
in fc |2.5| 

7 Newly formed stellar particles can move quickly between simu¬ 
lation outputs during the star formation burst, owing to feedback 
induced mass redistribution, so we integrate star formation within 
this larger radius. 

8 Results averaged over O.lGyr or shorter time-scales are qualita¬ 
tively similar but show more rapid fluctuations in sub-components 
of longer bursts so we selected 0.2Gyr for the sake of clarity. 
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Each strong burst of star formation reduces central den¬ 
sity of DM, so the final slope and core size are a consequence 
of several burst episodes over a Hubble time. The overall 
effect is small in mlO, because the star formation rate de¬ 
creases to very low values at 2 < 4. From the comparison of 
DM only and feedback simulations it is clear that a small 
difference in central DM concentration was established early 
on and stays largely unchanged until late times. 

A small, fluctuating and early decrease in the DM 
concentration is also seen in more massive halos, e.g. 
ml0hl297 and mil. However, the amount of dark mat¬ 
ter in a central region only decreases significantly once the 
central region finishes rapid growth, at z ~ 1 in ml0hl297 
and 2 ~ 3 in mil. After this stage, DM only simulations 
show an approximately constant amount of dark matter in 
the central region while baryonic simulations successfully 
evacuate a large amount of dark matter from the center. 
Unlike mlO, they have several ongoing bursts of star for¬ 
mation after the rapid-buildup stage, each of which removes 
a significant amount of dark matter. It appears that having 
strong bursts of star formation and outflows after the inner 
halo buildup slows down is the key to produce a long lasting 
shallow DM density profile. At early times, during the fast 
buildup stage, shallow profiles are not fully established as 
showed in Figure [5] 

In ml2i and ml2v, removal of DM after the peaks of 
star formation is also seen at 2 > 1 when fluctuations are 
large. After 2 ~ 1, the star formation continues at modest 
level without rapid bursts. At the same time, the enclosed 
DM mass grows slowly. In these massive halos re-accretion 
of dark matter in the center occurs when the star formation 
rate is low and hierarchical assembly is slow. To explain this 
we followed the central accumulation of baryonic material 
and found out that in both halos baryons start to dominate 
central mass at 2 ~ 1.5. As a consequence dark matter gets 
contracted, increasing the amount in the inner halo (see §[4|. 
This effect is stronger in the ml2i simulation that is more 
baryon dominated, and accumulates baryons faster at late 
times. |f] 

The core forms after multiple starbursts rather than 
one single blow-out, which is consistent with the mechanism 


discussed in Pontzen & Governato (2012 I (see also Ogiya & 
Mori|2014). After a blow-out some re-accretion of DM does 


occur but even after several Gyrs the amount of DM enclosed 
in central region does not return to the pre-burst level, indi¬ 
cating that the effect on DM is long lived. These trends are 
the most obvious at A/h ~ 10 11 Mg, i.e. in mil. However, 
in more massive systems we see that cuspier profiles are re¬ 
built at later times. While star formation proceeds to late 
times, it is often spread throughout the disk. In a companion 
paper (Muratov et al.i|2015l we show that at late times the 
star formation activity is not able to eject large quantities 
of material from galaxies, which is why this continuous star 
formation does not “heat” and remove DM from the center. 
Even if the assembly of central regions of the dark matter 


9 In ml2v the star formation rate is low and the amount of 
baryons in the central region changes very slowly in the final 
several Gyrs. The enclosed dark matter at those times is also 
affected by ongoing minor mergers, that cross very close to the 
center and are later disrupted causing variations in the enclosed 
density. 


halos slows down at relatively early times, increase in central 
concentration of baryons at late times can rebuilt DM cusps 
via adiabatic contraction. Our simulations suggest that the 
contribution of minor mergers to the re-growth of cusps (e.g. 
Dekel et al.|2003 Laporte fe Penarrubia|2014 1 at late times 


is likely sub-dominant or negligible, but a larger number of 
halos is needed to confirm this for halos with diverse growth 
histories. 

We have already seen that core formation depends 
strongly on halo mass but also on the presence of signifi¬ 
cant star formation episodes after the central region stops 
accreting dark matter. This is critical at halo masses where 
large cores start to develop. We showed this directly by us¬ 
ing the initial condition of mlO, simulated with a slightly 
different local feedback coupling scheme that increases late 
time star formation and its burstiness: the outcome is a 
larger DM density core and a shallower center slope com¬ 
pared to the original FIRE mlO simulation analyzed here 
(see Onorbe et al. 2015 I This means that if there are bursts 
of star formation activity occurring in dwarf galaxies/halos 
around this mass at late times they could result in shallow 
density profiles by present time. 

The star formation history of mlO is significantly dif¬ 
ferent from ml0hl297. It forms most of the stars before 
2 = 2 and becomes passive at late time whereas other halos 
remain actively star forming at present. From the obser¬ 
vations of dwarf galaxies in the Local Group, |Weisz et al.[ 
(2014) show that a large fraction of dwarf galaxies have ac¬ 
tive late time star formations, especially for field galaxies 
and galaxies with A/* > lO b M 0 . Therefore, the star for¬ 
mation histories of our simulated galaxies are well within 
the observed range. Our limited statistics suggest that at 
A/* ~ 10 6 — 10 7 Mq SFR history is closely linked to forma¬ 
tions of large cores (see also Onorbe et al.|20151. 


3.5 Halo Expansion or Baryonic Contraction? 

In the previous section, we investigated the shapes of pro¬ 
files of our simulated halos. Here we examine the effect of a 
central concentration of baryons on the dark matter profiles. 

First, we examine the net effect of halo expansion via 
feedback and halo contraction owing to central concentra¬ 
tion of baryons. Feedback is dominant in shaping flatter pro¬ 
files in lower mass halos, but baryonic contraction largely 
cancels the feedback effect in Milky Way mass halos, such 
that their final profiles are only slightly shallower than the 
NFW profile. 

In order to estimate the effect of baryonic contraction 
on the dark matter profiles, we follow Blumcn thal et ah] 
(1986) to calculate the final dark matter mass distribution 
M x , given the final baryon mass distribution mb and the 
initial total mass distribution Mf. 

r[m b (r) + M x (r)\ = riMi(n) = nM x (r)/( 1 - F b ), (3) 

where F b is the fraction of dissipational baryons, and r 
are initial and final orbital radii respectively. We follow the 
simple semi-analytic model (Dutton et al. 2007|, and assume 
that initial total halo density profile is the one from our DM 
only simulation, and use the final distribution of baryons 
(stars and gas) in the full physics simulation to estimate 
the contraction. Therefore, we use m b from hydrodynamical 
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Figure 5. (Left) Dark matter density slope, a, within l — 2%R vu -(z = 0), RvO, as a function of redshift/cosmic time in mlO, ml0hl297, 
mil, ml2v and ml2i respectively (from top to bottom). The slope is measured at fixed physical radius at all times. Blue lines shows 
the variations on the time scale of our simulation output (typically 20-30 Myrs) while the black lines show the average values over 
0.5Gyr. (Right) Time evolution of the enclosed dark matter mass within 0.02R v o f° r DM only (black dash-dotted) and hydrodynamical 
simulations with feedback (blue, dashed) and star formation rate within 0.1R v o (red solid), all averaged over 0.2Gyr. Green vertical lines 
in mil panels show the times of close passages of a subhalo. Cores form when the central DM accretion stops but the star formation is 
still bursty, as seen in ml0hl297 and mil. 
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simulations with stellar feedback, and M; is the halo mass 
from dark matter only simulations. □ 

Figure [ 5 ] shows density profiles from collisionless simu¬ 
lations, full feedback simulations, and from models with a 
contracted DM halo. Left panels show three different DM 
profiles for each halo and baryon density profiles from feed¬ 
back simulations for reference. Right panels show the total 
density profiles including DM and baryons (gas and stars). 

In mlO, the estimated effect of adiabatic contraction is 
very small because the fraction of baryons in the center of 
mlO is small. The baryon density is less than one tenth of 
dark matter density near the center. However, the feedback 
slightly expands the DM and forms a small core. 

In ml0hl297 and mil, stellar feedback strongly af¬ 
fects dark matter distribution, making it much shallower 
than in the corresponding collisionless run. It is interest¬ 
ing that in ml0hl297, feedback significantly flattens the 
DM profile, creating a large core, but the halo remains dark 
matter dominated at all radii. In mil, the baryon density 
is a slightly larger fraction of the total, but still significantly 
lower than the corresponding dark matter density in the 
collisionless simulation. The contracted profiles for both of 
these halos are therefore very similar to the dark matter 
profiles from the corresponding collisionless simulations. 

In more massive halos, the ml 2 series, baryonic con¬ 
traction is expected to significantly steepen the DM profiles, 
because their central regions are baryon dominated. While 
ml2v shows a strong effect of feedback, the profile of ml2i 
is relatively steep in the center, similar to the NFW profile. 
However, when compared to the expectations of baryonic 
contraction we see that the resulting profile is much shal¬ 
lower than contracted NFW halo for all of the plotted simu¬ 
lations. This demonstrates that even in these massive halos, 
feedback has a strong effect on shaping the final dark mat¬ 
ter profiles, and largely cancels out the effect of contraction. 
In general the expansion of the halo by the stellar feedback 
causes an order of magnitude difference in the density pro¬ 
files around 1 kpc in 10 11_12 Mg halos when compared to 
expectations of a simple baryonic contraction model. 

While feedback effects on the DM distribution are sub¬ 
stantial even in ml2 series halos, in those halos baryons 
completely dominate the central few kpc at z=0. This is 
why the differences between the total matter density pro¬ 
files in simulations with feedback and profiles expected from 
the contracted original DM halo are relatively small (right 
panels). We will return to these total matter profiles shortly 
and show that even this small effect has measurable con¬ 
sequences in the circular velocity curves of galaxies. It is 
interesting to note that the total matter distribution in the 
inner 20% of the R V1T of ml2 series is well approximated by 
the isothermal density profile (p oc 1/r 2 ). 

The results from our ml2v simulation are consistent 
with the strong feedback run in Maccio et al. (2012b I who 
showed core formation in a 7 x 10 ii Mg halo. We do how¬ 
ever find slightly higher central density of dark matter than 


10 We also tried to account for the loss of baryons due to feedback 
but the mass loss is much smaller than the total mass, especially in 
ml2 series. The difference in circular velocity with or without the 
missing baryons is less than a few percent in mil and negligible 
at higher masses. 


reported in Maccio et al. (2012b). While this might be just 
a matter of small number statistics, our simulation results 
should give more accurate predictions for the central profiles 
both because of more realistically implemented stellar feed¬ 
back and because of the higher resolution. At masses similar 
to the Milky Way (Afh ~ 10 12 Mq), the dark matter density 
distribution in the center is only slightly shallower that the 
NFW profile as strong feedback effects and adiabatic con¬ 
traction of such expanded halo nearly cancel out. 


4 OBSERVABLE CONSEQUENCES 

The following subsections show how our simulations with 
stellar feedback can alleviate the tensions between previous 
simulations and observations, including the “lack” of bary¬ 
onic contraction, the “Too Big To Fail” and cusp/core prob¬ 
lems. 


4.1 Rotation curves and the Tully-Fisher 
normalization 

The distribution of matter in galaxies can be measured with 
the rotation curves. For disk galaxies, there is a tight rela¬ 
tion between their luminosity (or mass) and their circular 
velocity, so called the Tully-Fisher < |TuIly fe Fisher 1976) re¬ 
lation (TFR.). Here we examine the effect of stellar feedback 
and baryonic contraction on rotation curves. 

Figure[7]shows the rotation curves of halos from simula¬ 
tions and baryonic contraction calculations for mil, ml2v 
and ml2i. Here V c = ^GM(r)/r. The rotation curves of 
simulated galaxies do not show profiles expected from the 
NFW halos affected by adiabatic contraction. Instead, galax¬ 
ies have lower mass concentrations in the centers, resulting 
in lower circular velocities. Therefore feedback effectively 
prevents buildup of high densities expected from strong adi¬ 
abatic contraction. As expected, the rotation curves in the 
m 12 series are approximately flat at radii larger than sev¬ 
eral kpc, while for the massive dwarf mil the rotation curve 
is rising. 

The effect of feedback, which cancels the effect of bary¬ 
onic contraction, turns out to be very important for the nor¬ 
malization of the TFR. We do not perform a detailed com¬ 
parison with observations, which would require mimicking 
observational measurements of the rotational velocity and 
luminosity. We plan to study this in future work. Here, in¬ 
stead, we focus on the relative effect of the feedback on dark 
matter profiles that determine the normalization of TFR. To 
show this effect and compare our simulated galaxies with the 
observed TFR, we measure the circular velocity of the halo 
V c = y/GM(r)/r. We measure circular velocity at 2.2 ‘disk 
scale length’ V 2 . 2 , approximately mimicking frequent obser¬ 


vational approach (Dutton et al. 2010). We first measure the 


half mass radius of the stellar distribution rr /2 and use the 
relation between half mass radius and scale radius of an ex¬ 
ponential disk to define disk scale length as ra = ri/2/1.67. 

Figure [8] shows the TFR of main galaxies with 10 9 < 
M,/Mq < 10 11 , and the best fit of the observed TFR from 
|Dutton et al.| (|2010[). They derived the TFR using the data 
from Courteau et al. (20071 and with the best fit: 
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Figure 6. (Left) Dark matter density profiles of mlO, ml0hl297, mil, ml2v and ml2i at z = 0. Different line colors show the 
expected DM profile in simulations with baryons and feedback (solid; red), in collisionless simulations (dashed; black) and in calculations 
including baryonic contraction (dot-dashed; blue). The green dotted line shows the total baryon density, including both gas and stars, 
in feedback runs. (Right) Total density profiles of the same halos, including both dark matter and baryons. The Power convergence radii 
are shown as dashed vertical lines. In halos where baryons dominate in the central regions, total and dark matter densities based on the 
simple baryonic contraction model are higher than the actual densities in our simulations with baryonic feedback. Feedback effectively 
cancels the effect of baryonic contraction. 
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This relation was derived for relatively massive galaxies, 
most with logV c [km/s] > 1.8. Similar relations were found 
for dwarfs, however with significantly enlarged scatter and 
non-uniform way of measuring V c (e.g. Ferrero et al.|20l2 l. 
We therefore limit discussion here to galaxies with M* > 
1O 9 M 0 . 

It is clear that the strong feedback which reduces the 
effect of adiabatic contraction is a necessary ingredient in re¬ 
producing and explaining the TFR in massive galaxies with 
M* > 10 1 o Mq. While there are direct effects of feedback 
on the distribution of baryons within galaxies feedback ef¬ 
fect on the distribution of dark matter is also an important 
ingredient in establishing the TFR. Simulated galaxies ap¬ 
pear to better match the observed TFR than our model with 
baryonic contraction. The circular velocities in baryonic con¬ 
traction calculations are higher by a factor of 1.2-1.5 than 
V c in simulations. 

Our findings confirm previous conclusions that the lack 
of effective contraction is necessary to explain the Tully- 


Fisher relation (Dutton et al. 2007 Maccio et al. 2012b). 


This also explains why previous generations of simulations 
without efficient feedback had trouble matching the normal¬ 
ization of TFR (e.g. Steinmetz & Navarro 2000). In these 
models too much gas collapsed to the center exerting strong 
contraction of DM halo without previously affecting the DM 
distribution. 


4.2 Implication for the ‘Too Big To Fail’ problem 


In addition to the cusp/core problem, cold dark matter 
simulations are also challenged by another problem, the so 


called “Too Big to Fail” problem (Boylan-Kolchin et al. 2011 
|Garrison -KimmeI et al. 2011}. The simplest version of this 
problem is that the observed Milky Way’s satellite galaxies 
have much lower central circular velocities than sub-halos 
from cosmological collisionless cold dark matter simulations. 
This either means that massive sub-lialos do not have corre¬ 
sponding match in observed satellite galaxies or that central 
regions of predicted cold dark matter sub-halos are too dense 
compared to observed halos. This seems to be generic prob¬ 
lem, independent of the halo formation history as similar 
effects are also observed in the Local Group and for dwarf 


galaxies in general (Ferrero et al. 

2012 Garrison-Kimmel 

et al. 2014 Papastergis et al. 2015 

) including non-satellite 


galaxies. 

We have already shown that stellar feedback can re¬ 
duce the central dark matter density, which could poten¬ 
tially resolve this discrepancy, without invoking different a 
type of dark matter. In Figure [9] we show the circular ve¬ 
locity profiles of the central kpc of mlO, ml0hl297 and 
ml0hll46 at z=0 along with their corresponding dark 
matter only simulations. These are our best resolved sys¬ 
tems with galaxy stellar mass ~ 2 x 10 6 — 8 x 10 7 Mq at 
z=0, which are close to the stellar masses of the galaxies 
for which the “Too Big to Fail” problem was demonstrated 
(Boylan-Kolchin et aJ7||201l| Garrison-Kimmel et al. 20131. 
For comparison we also include the observational data from 


Milky Way satellite galaxies (Strigari et al. 2007 Walker 
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Figure 7. Rotation curves of halos, including mil, ml2v and 
ml2i. Red solid lines represent rotation curves from simulation 
with feedback, black dashed lines show results from DM only sim¬ 
ulations, whereas blue dot-dashed lines represent rotation curves 
with baryonic contraction. Dashed lines show the region within 
the Power radius. Their halo masses are shown in the brackets. 
The rotation curves in simulations with baryonic feedback are 
lower than in a simple baryonic contraction model. The simu¬ 
lated Milky Way-mass halos show flat rotation curves, while the 
large dwarf galaxy mil shows a rising rotation curve. 


|et al. J2009 Wolf et al. J2010) and Local Group field galaxies 
( Kirby et al.|2014 l. It is clear that feedback strongly reduces 
circular velocities in the central few hundred pc with respect 
to collisionless cold matter simulations, in some cases by a 
large factor. Such reduced circular velocity implies that ob¬ 
served dwarf galaxies (including satellite galaxies) should 
not be associated with halos and sub-halos from DM only 
simulations with the same circular velocity, instead these 
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Figure 8. Stellar mass — Tully-Fisher relation for observed 
galaxies (Eq. and for the simulated halos (mllh383, mil, 
ml2v, ml2q and ml2i). The y-axis is the rotation velocity mea¬ 
sured at 2.2 disk scale length. The shaded regions show Eq. [~i] 
with one sigma uncertainty (a = 0.039) (cyan) and two sigma un¬ 
certainty on the zero-point (grey) respectively. Stellar feedback, 
which counteracts the effects of adiabatic contraction, appears 
necessary to establish the observed normalization of the Tully- 
Fisher relation. 


should be connected to the predicted higher circular veloc¬ 
ity analogs, whose circular velocity is now reduced owing to 
feedback. Our results strongly suggest that this effect would 
dramatically reduce the number of “massive failures” and 
can alleviate or potentially solve the “Too Big to Fail” prob¬ 
lem. Our findings qualitatively agree with hints in previous 
work ( Brook &; Di Cintio|2015| . 

It is interesting that the high stellar mass dwarf galaxies 
(e.g. ml0hll46 and ml0hl297) have more significant re¬ 
ductions in central rotation velocities (and thus dynamical 
masses) compared to low stellar mass dwarf galaxies (e.g. 
mlO). This causes the rank order of V c at small radius (e.g. 
500pc) not to correspond to the rank order of their V m ax 
or their stellar mass, as illustrated in the middle and lower 
panels of Figure [9] Direct comparison between dark mat¬ 
ter only simulations and simulations with baryons is even 
more complex, and rank order matching of V c , or V pea k from 
measurements at small radii might lead to incorrect physical 
interpretations. 

Our results only indirectly address the “Too Big To Fail” 
problem in satellite galaxies, because the galaxies we con¬ 
sider here are not satellites but field galaxies. Satellite galax¬ 
ies of relevant mass in our ml2 simulations do not have 
sufficient mass resolution to study their dark matter distri¬ 
butions at the galaxy centers. The effect of host galaxies on 
satellites, e.g. tidal stripping, could also modify the structure 
of DM halos ( Zolotov et al.|[2012 Brooks fe Zolotov||2014 |. 
However, the effect of lowering the circular velocity of galax¬ 
ies is generic in the range of stellar masses 2 x 10 6 — 3x 10 9 Mq 
and we therefore believe that satellite galaxies in this mass 
range will be affected in the same systematic way. 


4.3 Central Slopes and Core sizes 

pT et al. (2015} and Oh et al. (2015) measured detailed 
density profiles of field dwarf galaxies, and found flat DM 




Figure 9. Rotation curves illustrating the TBTF problem, plot¬ 
ted over a range of radial scales. We have included halos mlO 
(M* = 2.3 X 10 6 M o; black, thin), ml0hl297 (M* = 1.7 X 
10 7 Mq; blue, normal) and ml0hll46 (M* = 7.9 X 1 O 7 M0; 
green, thick) and their corresponding DM only simulations. Thick 
dashed lines represent the halos from the collisionless simulations 
while thick solid lines represent the halos from the hydrodynam- 
ical simulations with feedback. The thin lines show the velocity 
curves at radius smaller than the Power convergence radius. Black 
squares show the data from Milky Way bright satellite galaxies 
( [Strigari et al.|2007| | Walker et al.|2009| |Wolf et al.|2010| ) , while 
red circles show the isolated dwarf galaxies in the Local Group 
( [Kirby et al.|2014| . The panels show three different scales of the 
same plot to illustrate that the order of rotation curves at small 
scales may not imply the order of halo masses. Stellar feedback 
reduces the central circular velocity such that the rotation curves 
can match the observed dwarfs, suggesting that baryonic feedback 
may solve the “Too Big To Fail” problem. Observational errors are 
typically smaller than a few km/s (not shown for clarity). 
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Figure 10. Slope of dark matter density profile, a, from our 
simulations (measured at r = 0.3 —0.7kpc) compared with a from 
observations (typically measured at a few hundred pc). Blue solid 
circles represent the simulated halos at z = 0. Red hollow squares 
represent the observed dwarf galaxies from THINGS ( |Oh et ah] 
1 2011 1 [Walter et al.|2008] |, whereas black hollow triangles represent 
the dwarf galaxies from LITTLE THINGS < |Oh et al.|2015|[Hunter| 
et al. 2012 j ). In the overlapping mass range, the simulated dwarf 
galaxies have central DM profile slopes in good agreement with 
the observed dwarfs. 


profiles near the center, in contrast with the cuspy NFW pro¬ 
files, predicted from N-body simulations. In massive galax¬ 
ies, such as the Milky Way, baryons dominate in the central 
few kpc (e.g. |Courteau fe D utton| j2015|l making measure¬ 
ments of central dark matter properties extremely difficult. 
We therefore focus on lower mass galaxies/halos with reli¬ 
ably measured central DM profiles. 

In Figure |10[ we show that FIRE halos at z = 0 are 
in good agreement with the observed slopes of central DM 
density profiles (see also |Governato et al7||2()l 2[ ) . This sug¬ 
gests that the inclusion of stellar feedback in our simula¬ 
tions helps resolve so called “cusp/core” problem observed 
in low mass galaxies. However, the observed scatter is large 
and number of simulated objects is limited in the observed 
mass range. It is therefore clear that much larger sample of 
model galaxies/halos as well as more detailed accounting for 
methodology and selection used in observations are needed 
to test our model in detail. 

Figure m shows the DM core radii of our simulated 
galaxies as a function of their baryonic masses and compares 
them to observations from Walter et al. (2008). In general, 
for m b ~ 10 7 — 10 10 Mq E3 , core radius increases with bary¬ 
onic mass. There is a broad agreement with observed core 
sizes. A more detailed comparison will require a larger num¬ 
ber of halos and more detailed modeling of the methodology 
used to calculate the core sizes in observed galaxies. 


11 m is defined as the total baryonic mass within 20% /i. v i, of 
the halo. 
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Figure 11. Comparison between the core radii of THINGS dwarf 
galaxies and our simulations as a function of baryonic mass 
(see footnote \n) at z = 0. Red hollow squares show the observed 
core radii | |Walter et al.|[2008| |Oh et aj~]|20 11 | i while blue circles 
show the core radii from our simulations. Within the plot range, 
the core size increases with baryonic mass, which is largely con¬ 
sistent with the observed cores. A larger number of observed and 
simulated dwarfs is necessary to draw stronger conclusions. 


5 DISCUSSION 

We use the FIRE suite of high-resolution cosmological 
“zoom-in” simulations with physical feedback models to 
study the properties of the central regions of dark matter 
halos. Our simulations have higher resolution (in most cases 
the highest resolution to date at a given halo mass), and 
more explicit and comprehensive implementation of stellar 
feedback than simulations used previously to study bary¬ 
onic effects on DM profiles. Critically our models contain no 
“freely adjusted” parameters tuned to any particular results. 
We characterize the evolution of halo properties and corre¬ 
lation with galaxy and halo mass and explore the effects of 
star formation driven feedback and adiabatic contraction on 
the slopes, cores size and circular velocity profiles of galactic 
halos. 

We find that at 2 = 0 the central slopes of DM den¬ 
sity profiles measured at ~ l%_R v ir are shallow in ha¬ 
los Mh ~ 10 1() — 10 11 Mq, but the slopes are cuspier at 
lower and higher masses. We see a sharp transition around 
Mh ~ 10 10 Mq from cuspier profiles at lower masses to shal¬ 
low, core-like profiles at higher masses. Efficient feedback 
continues to at least Mh ~ 10 1j Mq where the core-like pro¬ 
files that form during earlier evolution are contracted by 
baryons at late times into steeper profiles. Final profiles are 
similar or flatter than NFW, and therefore significantly flat¬ 
ter than expectations for a contracted NFW halo. 


Our results are in broad agreement with others found 
in the literature. For example we find that feedback is ef¬ 
ficient in forming large cores at halo masses ~ 10 10 Mp)- 


few x10 11 Mq which is similar to previous findings ( 

Gov- 

ernato et al. 2012) Di Cintio et al. 2014; |Pontzen & Gov- 

ernato 

20141. It appears that simulations with bursty star 


formation and outflows, but different small scale feedback 
implementations, affect dark matter profiles in a qualita¬ 
tively similar way. However we also found some interesting 
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A U (l%core) 
(erg) 

AU (3% core) 
(erg) 

E sn 

(erg) 

m09 

5.85e52 

6.45e53 

2.68e53 

mlO 

1.08e53 

2.02e54 

2.13e55 

mil 

3.43e55 

6.46e56 

2.22e58 


Table 3. Comparison of the difference in gravitational potential 
energy AU of halos from dark matter only simulation and the 
same halos with constant density core of 1% of _R v ir (left) or 3% 
of -R v ir (middle), and their total supernova energy (E S n, right), 
obtained from the corresponding simulations with feedback. 
Energy input from supernovae alone is sufficient to produce 
cores in mlO and mil. 


differences. This is not surprising, given that our simulations 
are in most cases higher resolution and have explicit treat¬ 
ment of feedback on small scales. For example we find that 
in the halos with large cores, cores are established at around 
z ~ 1 and can grow even at low redshift while other authors 
find that cores formed very early already at z = 2 — 3 (e.g. 
Madau et al.||[2014 ). This could be a consequence of differ¬ 


ent star formation histories that can change core formation 
| |Onorbe et al.| |2015). For example, these authors note that 
earlier simulations or dwarf galaxies with “sub-grid” feed¬ 
back models produced too many stars at a given halo mass. 
Alternatively, these differences could indicate that different 
treatments of material affected by feedback (cooling preven¬ 
tion vs. explicit model) cause differences in DM profiles. 

We also find cuspier DM p rofiles at around 6 — 7 x 
10 11 Mg than the one reported in Maccio et al. (2012b) and 
therefore weaker halo expansion and somewhat steeper de¬ 
pendence of slope on M*/Mh than the relation presented 
Di Cintio et al. (2014). Overall, better statistics from a 


larger number of simulated halos are needed for a more ro¬ 
bust analysis of these differences. 


5.1 Energetics 


Given that the efficiency of conversion of halo gas into stars 
increases from dwarf galaxies to massive galaxies, it is natu¬ 
ral to connect feedback effects to the energy available from 
stellar feedback. We therefore compare the amount of energy 
available from feedback with the energy needed to overcome 
the gravitational potential and move dark matter outside 
of the central region. The simplest estimate is to calculate 
how much energy gets injected into the ISM from the SNe 
only, which represents a lower limit to the available energy 
budget. 

Previously, |Gnedin &; Zhao|(|2002| ), [0giya fe Mori| ( [2C)ll| ) 
and Garrison-Kimmel et al. (20131 claimed that the to¬ 


tal feedback energy released in SNe is insufficient to re¬ 
move enough dark matter to form large cores. In particu¬ 


lar, Garrison-Kimmel et al. (2013) tested if stellar feedback 


can lower central densities to a degree needed to explain 
the “Too Big To Fail” problem. They simulated supernova 
feedback with time-varying potential and found the number 
of SNe needed to match observed profile of a halo hosting 
M* ~ 10 6 Mq galaxy exceeds the number of SNe produced in 
most of the dwarf galaxies for the typical initial mass func¬ 
tion. However, they did not consider the full growth history 
of the halo. The frequent mergers and star formation bursts 
at higher redshift could dynamically heat up more dark mat- 
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ter in the center, when halo was smaller. Furthermore, they 
assume mass loading of SNe driven winds smaller than we 


find for FIRE galaxies of comparable masses (Muratov et al. 


20151. Finally their selected halo has a high concentration 


compared to a typical sub-halo that is expected to host ob¬ 
served dwarf galaxies. As we show below, there is a sufficient 
amount of energy available to couple to dark matter in the 
relevant halo and galaxy mass regime. 

Table [3] shows the energy needed to create a constant 
density core with radius 1% or 3%_R v ir in halos from our 
dark matter only simulations and the total supernova en¬ 
ergy inferred from our hydrodynamic simulations by z=0. 
We constructed a constant density core in our “cuspy” DM 
only simulations by keeping the density profile outside the 
core unchanged and moving the excess mass within the core 
radius to infinity. Then we calculate the total potential en¬ 
ergy for the initial and cored profiles with the formula below 


U = -4?rG 


r 


M(r)p(r)r&r, 


(5) 


and report the difference, A U, in Table [3] 

To estimate the supernova energy, we assume the energy 
from one supernova E sn = 10 51 erg, the fraction of massive 
stars which can produce supernovae £(m* > 8Mq) = 0.0037 
(Kroupa 20021, and the mean stellar mass M mean = O.4M 0 . 
The total supernova energy of a halo is given by 


-Esn = M„/(< m* >) * £(m* > 8M 0 ) * E SI 


( 6 ) 


From Table [3] we see that for m09, the amount of SNe 
energy is not sufficient to create a large core at 3% R v n even 
if it all couples (via secondary gravitational interactions) to 
the DM. Even creating 1% R v i r core is difficult as it requires 
that more than 20% of the available energy is coupled to 
dark matter, which is unlikely given the indirect connection 
via the change of gravitational potential and that a large 
fraction of this energy is in heavily mass loaded winds that 
move rapidly ( Muratov et al.| 2015l. 

However, for mlO, a small core is energetically possible 
and even a 3%7? v ir core requires less than 10% of the avail¬ 
able energy. We see some signs of profile flattening within 
the inner 1% of R v i r , but not a fully developed core. How¬ 
ever, for the same halo mass Onorbe et al. (20151 show that 
a slightly different star formation history can cause a much 
larger effect and form a central core with a radius of 1-2 

%-Rvir- 


In mil, the supernova energy is three orders of magni¬ 
tude higher than what is needed to create a small core at 1 
% of Rvir and even 3% of R v i r core can be created with few 
percents of coupling efficiency. It is therefore not surprising 
that this halo indeed hosts a large core. Even though the 
depth of the potential well increases in more massive halos, 
the amount of stellar mass is a steep function of halo mass, 
which provides sufficient energy to affect central dark mat¬ 
ter profiles. This is why we see a relatively sharp transition 
at ~ 10 10 M Q from low mass cuspy halos to core-like halos 
at higher masses. 


5.2 Correlation with star formation 

Here we show a simple connection between bursts of star 
formation that cause strong gas ejection episodes, and the 
change in the amount of dark matter in the central part of 
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halos. In Figure [T2} we plot the change in the enclosed dark 
matter mass as a function of the peak star formation rate 
in mil, the halo in which the effect of feedback on dark 
matter is one of the strongest in our sample. We focus on 
bursts with high peaks of star formation and neglect low 
star formation rate episodes as they typically do not show 
strong outflows (Muratov et al. 20151. In Figure 12 we plot 
the rate of change in the central amount of dark matter as a 
function of the peak star formation rate. The star formation 
rate and mass change are measured as averages over 0.2 Gyr. 
The rate of change of the DM mass is 


AM _ radm(t + 0.2Gyr) - m dm (t) 
At W “ 0.2Gyr 


(7) 


where mdm(t) is the enclosed DM mass within 2%7? v o (the 
virial radius at 2 = 0) averaged over 0.2Gyr and t is the 
beginning time of the interval. We measure star formation 
between t — 0.2Gyr and t, within a larger radius of 0.1 R v o to 
make sure to include stars that move out of the very center 
between the two time intervals. [3 

A significant decrease of mass enclosed in the central 
2%R v o occurs just after a strong starburst when a large 
amount of dark matter is “heated” and effectively pushed 
out of the central region. At the corresponding times DM 
only simulation does not show negative mass change, except 
in the case of the largest burst which is triggered by the 
close passage during a merger (as indicated in Figure[5|, i.e. 
the merger dynamically alters the profile. Hence we conclude 
that there is a correlation between mass removal from the 
center and star formation rate. 

This correlation suggests that a strong burst of star 
formation can provide sufficient feedback to remove a sig¬ 
nificant amount of baryons, which then cause a decrease in 
the central potential and lower the concentration of DM. 
This scenario is consistent with the mechanism suggested in 


Pontzen & Governato (2012) in which repeated changes of 


central gravitational potential transfer energy to the orbits 
of DM particles causing a central density decrease. 

Some fraction of the removed DM does return, i.e. cores 
get partially rebuilt, however when large cores are present 
this is a relatively small effect. When cores are established 
after the period of rapid early halo growth, they survive 
largely intact for at least several Gyrs, before a new major 


bursts of star formation and gas expulsion sets in (see § 3.41. 


In more massive halos (ml2), shallow profiles that form at 
higher redshift steepen by redshift zero significantly. In these 
halos we see continuous late time star formation without sig¬ 
nificant galactic wind episodes (Muratov et al. 2015) which 
continuously increases central density of baryons that dom¬ 
inate the central potential. As a consequence dark matter 
halos are pulled inward by adiabatic contraction changing 
shallow profiles into cusps. 


12 We have also tried to average SFR over 0.1 Gyr and considered 
the delay of 0.1 Gyr. The enclosed DM mass drops 0.1 Gyr after 
the peak SFR in away similar to Figure 12 (i.e. it is negative), but 
with the larger scatter and greater changes in DM mass. We note 
that this simple measurement cannot use much longer intervals 
for averaging of SF and for time delay because on large time scales 
the SFR bursts would be “washed out” and we could have multiple 
burst episodes within a long delay time interval. 
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Figure 12. The rate of change of enclosed dark matter mass 
within 2% of R v q as a function of peak star formation rate within 
10%-RvO hi mil (red circles) and dmll (black hollow squares) 
simulations, both averaged over 0.2 Gyr. In both cases we plot 
the SFR from the corresponding feedback simulations at the same 
cosmic time. The star formation rate is measured 0.2 Gyr ahead 
of the rate of change of enclosed dark matter mass. Strong bursts 
of star formation are followed by a reduction in the enclosed DM 
mass in the simulation with feedback. 


While a burst-driven core formation mechanism is con¬ 
sistent with our findings, we cannot exclude a contribution of 
other dynamical mechanisms, such as the motion of dense 
baryonic clumps within galaxies with respect to the halo 
centers. On the other hand, core formation via enhanced dy¬ 
namical friction from the dense infalling sub-structure (e.g. 
El-ZanFet al.|2001) is unlikely to play a significant role for 


the halo mass range probed here, because feedback lowers 
the density of infalling sub-lialos relative to their DM-only 
counterparts (see Figure [6] all of the infalling sub-halos in 
our sample have Mh < 10 M@). 


5.3 Significance for the dark matter detection 

The dark matter profile in the Milk Way is important in 
studies of indirect detection of dark matter particles from 
annihilation or decay signals. Recently, extended emission 
in gamma rays from the galactic center has been reported 
based on the data from the Large Area Telescope aboard the 
Fermi Gamma Ray Space Telescope (e.g. 

[enough 2011| ). To interpret this signal as 
annihilation of dark matter particles or to constrain its con¬ 
tribution it would be extremely useful to know the central 
dark matter profile (Abazaji an et akj |2014). The signal is 
consistent with the annihilation rate of thermally produced 


weakly interacting massive particle dark matter 

Goode- 

nough & Hooper|2009[|Hooper & Goodenough 2011 

Hooper 

& Linden 201l[Boyarsky et al. 2011; Abazajian & 

Xapling- 


Hooper & Good- 
a consequence of 
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hat||2012| |Gordon fe Macias 2013[ |Macias fc Gordon||2014[ 


Abazajian et al. 20141 although other possibilities remain 
valid alternatives (Abazajian 20111. 

While we have a small number of objects at the relevant 
mass, one robust finding from our simulations is that cen¬ 
tral density of cold dark matter will not correspond to con¬ 
tracted NFW profiles in Milky Way-mass halos. This helps 
constrain the range of values used in the modeling of the 
observed signatures. We find values of a ~ —1 to —1.4 at 
1 —2%7?vi r which is consistent with the best fit in |Abazajian| 
et al. (2014). However our results from ml2 halos suggests 


that deeper in the halo, the DM profile is likely shallower 
a ~ —0.7 to — 1.1 at < l%7? v j r . Given the slight differ¬ 
ences of the central slope definition and fitting procedures 
it would be interesting to test if the profiles shown here pro¬ 
vide a good match for the observed signal. We defer such 
more detailed comparison for future work. 


5.4 Limitations 


There are some clear limitations in our study of dark mat¬ 
ter halo properties. The number of simulated halos in our 
sample is limited because such high-resolution simulations 
are very time-consuming. In the future we plan to simulate 
a much larger number of halos to be able to extract direct 
statistical predictions and to compare to observed values, 
including the scatter in observations and theoretical predic¬ 
tions. 

Further limitation comes from the finite resolution. Our 
simulations are amongst the highest resolution cosmological 
simulations at z=0 to date and our m09 and mlO simula¬ 
tions can robustly resolve dark matter profiles on < 200 pc 
scales. However detailed comparison to observations requires 
converged results and the ability to exactly integrate circu¬ 
lar orbits over a Hubble time without N-body effects within 


300 pc in M h ~ 1O 1O_11 M 0 (M, ~ 1O 7 " 9 M 0 ) halos (Wal- 


|ter et ah]|2008[ ). In our ml0hl297, ml0hll46 and mil 

simulations, this is larger than many gravitational softening 
lengths, however the more stringent Power convergence cri¬ 
terion j ]2.4| ( |Power et al.| [2003) requires a further factor of 
several improvement in mass resolution to reach this limit. 

Resolution requirements can be even more dramatic if 
one wants to directly study the inner density profiles of small 
dwarfs that form within sub-halos of large spiral galaxies 
(e.g. to directly address TBTF problem). Such study would 
require simulating Mh = 1O 12 M 0 halos with > 10 9 parti¬ 
cles within the virial radius. We have shown, however, that 
many of the CDM “problems” can be explained by baryonic 
feedback in isolated, well resolved, dwarf galaxies suggesting 
a similar solution for satellite galaxies. 

We caution that the Power convergence criterion was 
derived for different time-stepping algorithm, different force 
accuracy and different softening values. Based on the DM 
profiles in Figure [T] we see that steep central slope in DM 
only simulations continue to at least factor of ~ 2 smaller 
radii than what is estimated by Power convergence radius. 
More importantly, direct resolution tests confirm the conver¬ 
gence in DM cusp profiles down to a factor ~ 2 — 3 smaller 
radii for our fiducial resolution. This suggests that this con¬ 
vergence criterion might be too conservative for our simula¬ 
tion setup. 

We did not study m!3 halo (Mh ~ lO lli M 0 ) from the 
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FIRE project sample because the star formation rate of 
this halo, at low redshift, is higher than the observed rates 
in observed galaxies of the same mass (H14). Some addi¬ 
tional physics, e.g. active galactic nucleus (AGN) feedback, 
is needed to explain the observations. Yet, the investigation 
of the profiles of these high mass halos is highly intrigu¬ 
ing both observationally and theoretically. Observationally, 
gravitational lensing provides an accurate measurement of 
the enclosed mass of those halos (Bolton et al. 20081, though 
dark matter only constitutes a minor fraction of the mass. 

Nevertheless, stellar kinematics and strong lensing do 


suggest cores in galaxy clusters with a ~ —0.5 (Newman 
et al. 2013b[a|. It is not clear if a single mechanism can ex¬ 


plain such shallow slopes over a large range. Mechanisms 
ranging from more frequent major mergers as well as pro¬ 
cesses such as AGN feedback, magnetic fields, anisotropic 
conduction, and cosmic rays, that are not yet incorporated 
in our simulations may be important in regulating late time 
star formation and affecting the core formation in massive 
halos (e.g. Peirani et al.|[2008 Martizzi et al.|2013 l. 

AGN feedback in particular might even affect halos with 
masses similar to our ml2 simulations, e.g. [Velliscig et al.| 
(2014) showed that AGN feedback can have a non-negligible 
impact on the halo properties (i.e. mass and profile) down 
to Mh ~ 5 x 1O 11 M 0 . These results suggest that the effect of 
AGN feedback, in addition to stellar feedback, could further 
lower the central density of the most massive halos in our 
sample. We note that regardless of the dominant feedback 
mechanism, the overall efficiency of feedback must be sim¬ 
ilar to what is seen in our simulations, as this efficiency is 
constrained by the observed M» — Mh relation. Our simula¬ 
tions provide a clean test for the effects of stellar feedback 
alone on the DM distribution. 


5.5 Dependence on star formation history 

|Onorbe et al. (2015]) compared our mlO simulation from 
H14 with the one with a slightly different supernova feedback 
coupling at smaller scales. In our default case, energy depo¬ 
sition is volume weighted while in the other version it was 
mass weighted. This creates slight differences in the feed¬ 
back and changes late time star formation. A 1-kpc core 
was formed in a halo with more prominent late time star 
formation, while our default mlO simulation shows a core 
at < 400pc and much higher central density of dark matter. 
It is likely that this strong sensitivity is caused by this halo 
mass being at the transition region from smaller to larger 
cores. 

This then suggests that when a large sample of sim¬ 
ulated halos is available, comparison to the observations 
around this mass scale will potentially help distinguish be¬ 
tween feedback models by analyzing their star formation 
histories and properties of their dark matter halos (Fitts et 
al., in preparation). For slightly more massive halos we see 
core formations in all cases we explored, regardless of their 
detailed SF history and the details of their feedback imple¬ 
mentations. 
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6 CONCLUSIONS 


7 ACKNOWLEDGEMETS 


We have explored cold dark matter profiles in simulations 
with stellar feedback. We used the FIRE suite of hydrody¬ 
namic simulations initially discussed in H14 and supplement 
these with 4 new dwarf galaxy simulations. We also run col¬ 
lisionless counterparts for all of these simulations. We show 
that baryonic simulations can successfully produce results 
consistent with observations and alleviate or solve several so 
called “problems” of the CDM: “cusp-core”, the lack of adi¬ 
abatic contraction or “halo expansion” and the “Too Big To 
Fail”, without any fine turning or introduction of adjustable 
parameters. Our main results are: 


1. The baryons have little influence on halos with Mh <S 
10 10 Mq because only a small fraction of available baryons 
are converted to stars, owing to feedback and the UV back¬ 
ground that suppress their star formation after the reioniza¬ 
tion ( Onorbe et al.|2015| Wheeler et al.|2015 l. The smallest 
halos are therefore perfect places for testing various theories 
of dark matter. 


2. The central slopes of dark matter density profiles are 
governed by halo mass and stellar mass. Profiles are shallow 
with relatively large cores for Mh ~ 10 10 — few xlO n M s 
and M* ~ 10 7 — few x10 9 Mq, where a ~ —0.5 — 0 and 
cores are r cor e > 1 kpc. Small central cores can also form at 
slightly lower masses M* ~ 10 6 Mq. This result is consistent 
with the observations of dwarf galaxies and can explain the 
“cusp-core" problem. 

3. Bursts of star formation and feedback start forming 
cores at early times but the cores are established typically at 
later times, e.g. in our mil simulation core is still growing 
at z < 1. Stable cores are established once central regions 
of halos stop their rapid growth. After this time (z < 2) 
removal of mass from the central region leaves a long term 
effect on the halo profile. We show that strong bursts of star 
formation are correlated with dark matter expansion. 

4. The total supernova energy in halos with Mh > 
1 O 1 O M 0 is sufficient to produce a core with radius 3%R v h, 
but not sufficient to make large cores in lower mass halos. In 
practice only a few percent of the available energy is trans¬ 
ferred into evacuation of dark matter from the central region. 

5. Baryonic contraction of dark matter halos becomes 
significant when central regions of halos are clearly dom¬ 
inated by baryons, which in our simulations occurs for 
Mh > 5 x 1 O 11 M 0 . However feedback in the progenitors of 
these massive galaxies significantly lowered central DM den¬ 
sity at z ~ 1 — 1.5. The cumulative effect of feedback and 
contraction is then a profile that in our ml2 runs is slightly 
shallower than NFW. This explains why the normalization 
of the Tully-Fisher relation requires no contraction or halo 
expansion with respect to the collisionless NFW profile. 

6 . Stellar feedback in galaxies with M* ~ few xlO 6 — 
10 s M© lowers the central density of DM when compared 
to dark matter only simulations and significantly reduces 
the rotational velocity near the center. This means that rel¬ 
atively low circular velocities of observed galaxies should 
correspond to much higher maximum circular velocities or 
virial velocities of halos or sub-halos in collisionless CDM 
simulation. This can solve or at least substantially alleviate 
the “Too Big To Fail" problem. 


We would like to thank N. Murray, J. Bullock, M. Boylan- 
Kolchin and P. Salucci for useful discussion. DK and TKC 
were supported in part by NSF grant AST-1412153, Heilman 
Fellowship and funds from the University of California San 
Diego. Support for PFH was provided by an Alfred P. Sloan 
Research Fellowship, NASA ATP Grant NNX14AH35G, 
and NSF Collaborative Research Grant ^1411920. EQ was 
supported in part by NASA ATP grant 12-APT12-0183, 
a Simons Investigator award from the Simons Founda¬ 
tion, and the David and Lucile Packard Foundation. CAFG 
was supported by NSF through grant AST-1412836, by 
NASA through grant NNX15AB22G, and by Northwest¬ 
ern University funds. The simulation presented here used 
computational resources granted by the Extreme Science 
and Engineering Discovery Environment (XSEDE), which 
is supported by National Science Foundation grant number 
OCI-1053575, specifically allocations TG-AST120025 (PI 
Keres), TG-AST130039 (PI Hopkins), TG-AST1140023 (PI 
Faucher-Giguere). 


REFERENCES 

Abazajian K., 2006, Pliys. Rev. D 73, 063513 
Abazajian K. N., 2011, JCAP, 3, 10 

Abazajian K. N., Kaplinghat M., 2012, Phys. Rev. D t 86 , 
083511 

Abazajian K. N., Canac N., Horiuchi S., Kaplinghat M., 
2014, preprint, (arXiv: 1402.4090) 

Agertz O., et al., 2007, MNRAS 380, 963 

Agertz O., Teyssier R., Moore B., 2009, MNRAS, 397, L64 

Barnes J. E., 2012, MNRAS, 425, 1104 

Begeman K. G., 1987, 

Blumenthal G. R., Faber S. M., Flores R., Primack J. R.., 

1986, ApJ, 301, 27 

Bolton A. S., Buries S., Koopmans L. V. E., Treu T., 
Gavazzi R., Moustakas L. A., Wayth R., Schlegel D. J., 
2008, ApJ, 682, 964 

Boyarsky A., Malyshev D., Ruchayskiy O., 2011, Physics 
Letters B, 705, 165 

Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, 
MNRAS, 415, L40 

Broeils A. H., 1992, PhD thesis, PhD thesis, Univ. Gronin¬ 
gen, (1992) 

Brook C. B., Di Cintio A., 2015, MNRAS, 450, 3920 
Brook C. B., Stinson G., Gibson B. K., Roskar R., Wadsley 
J., Quinn T., 2012, MNRAS, 419, 771 
Brooks A. M., Zolotov A., 2014, ApJ, 786, 87 
Bryan G. L., Norman M. L., 1998, ApJ, 495, 80 
Burkert A., 2000, ApJL, 534, L143 

Carlson E. D., Machacek M. E., Hall L. J., 1992, ApJ, 398, 
43 

Courteau S., Dutton A. A., 2015, ApJL, 801, L20 
Courteau S., Dutton A. A., van den Bosch F. C., 
MacArthur L. A., Dekel A., McIntosh D. H., Dale D. A., 
2007, ApJ, 671, 203 

Cullen L., Dehnen W., 2010, MNRAS, 408, 669 
Dave R.., Spergel D. N., Steinhardt, P. J., Wandelt, B. D., 
2001, ApJ, 547, 574 

Dekel A., Devor J., Hetzroni G., 2003, MNRAS, 341, 326 


© 0000 RAS, MNRAS 000, 000-000 




The Impact of Baryonic Physics on the Structure of Dark Matter Halos 21 


Del Popolo A., 2009, ApJ, 698, 2093 
Di Cintio A., Brook C. B., Maccio A. V., Stinson G. S., 
Knebe A., Dutton A. A., Wadsley J., 2014, MNRAS 437, 
415 

Dunstan R. M., Abazajian K. N., Polisensky E., R.icotti 
M., 2011, preprint, (arXiv: 1109.6291) 

Durier F., Dalla Vecchia C., 2012, MNRAS, 419, 465 
Dutton A. A., Maccio A. V., 2014, MNRAS, 441, 3359 
Dutton A. A., van den Bosch F. C., Dekel A., Courteau S., 
2007, ApJ, 654, 27 

Dutton A. A., Conroy C., van den Bosch F. C., Prada F., 
More S., 2010, MNRAS, 407, 2 
El-Zant A., Slilosman I., Hoffman Y., 2001, ApJ, 560, 636 
Elbert O. D., Bullock J. S., Garrison-Kimmel S., Rocha M., 
Oiiorbe J., Peter A. H. G., 2015, MNRAS, 453, 29 
Faucher-Giguere C.-A., Keres D., 2011, MNRAS, 412, L118 
Faucher-Giguere C.-A., Lidz A., Zaldarriaga M., Hernquist 
L., 2009, ApJ, 703, 1416 

Faucher-Giguere C., Keres D., Dijkstra M., Hernquist L., 
Zaldarriaga M., 2010, ApJ, 725, 633 
Faucher-Giguere C.-A., Hopkins P. F., Keres D., Mura¬ 
tov A. L., Quataert E., Murray N., 2014, preprint, 
(arXiv:1409.1919) 

Ferrero I., Abadi M. G., Navarro J. F., Sales L. V., 
Gurovich S., 2012, MNRAS, 425, 2817 
Flores R. A., Primack J. R., 1994, ApJL, 427, LI 
Garrison-Kimmel S., Rocha M., Boylan-Kolchin M., Bul¬ 
lock J. S., Lally J., 2013, MNRAS 433, 3539 
Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., 
Kirby E. N., 2014, MNRAS, 444, 222 
Gentile G., Salucci P., Klein U., Vergani D., Kalberla P., 
2004, MNRAS, 351, 903 

Gnedin O. Y., Ostriker J. P., 2001, ApJ, 561, 61 
Gnedin O. Y., Zhao H., 2002, MNRAS, 333, 299 
Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 
2004, ApJ, 616, 16 

Goodenouglr L., Hooper D., 2009, preprint, 

(arXiv:0910.2998) 

Gordon C., Macias O., 2013, Phys. Rev. D, 88, 083521 
Governato F., et al., 2010, Nature, 463, 203 
Governato F., et al., 2012, MNRAS, 422, 1231 
Hahn O., Abel T., 2011, MNRAS, 415, 2101 
Hooper D., Goodenough L., 2011, Physics Letters B, 697, 
412 

Hooper D., Linden T., 2011, Phys. Rev. D, 84, 123005 
Hopkins P. F., 2013, MNRAS, 428, 2840 
Hopkins P. F., 2014, ArXiv e-prints: 1409.7395, 

Hopkins P. F., Keres D., Oiiorbe J., Faucher-Giguere C.- 
A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 
445, 581 

Hunter D. A., et al., 2012, AJ, 144, 134 
Kaplinghat M., Keeley R. E., Linden T., Yu H.-B., 2014, 
Physical Review Letters, 113, 021302 
Katz N., White S. D. M., 1993, ApJ, 412, 455 
Kaufmann T., Mayer L., Wadsley J., Stadel J., Moore B., 
2006, MNRAS, 370, 1612 
Keres D., Hernquist L., 2009, ApJL, 700, LI 
Keres D., Katz N., Weinberg D. H., Dave R., 2005, MN¬ 
RAS, 363, 2 

Keres D., Vogelsberger M., Sijacki D., Springel V., Hern¬ 
quist L., 2012, MNRAS, 425, 2027 
Kim J.-h., Abel T., Agertz O., et al., 2014, ApJS, 210, 14 


Kirby E. N., Bullock J. S., Boylan-Kolchin M., Kaplinghat 

M. , Cohen J. G., 2014, MNRAS, 439, 1015 

Knebe A., Green A., Binney J., 2001, MNRAS, 325, 845 
Knollmann S. R., Knebe A., 2009, ApJS, 182, 608 
Kochanek C. S., White M., 2000, ApJ, 543, 514 
Kroupa P., 2002, Science, 295, 82 
Krumholz M. R., Gnedin N. Y., 2011, ApJ, 729, 36 
Laporte C. F. P., Penarrubia J., 2014, preprint, 
(arXiv:1409.3848) 

Leitherer C., et al., 1999, ApJS, 123, 3 
Loeb A., Weiner N., 2011, Physical Review Letters, 106, 
171302 

Lovell M. R., et al., 2012, MNRAS, 420, 2318 
Maccio A. V., Paduroiu S., Anderlialden D., Schneider A., 
Moore B., 2012a, MNRAS, 424, 1105 
Maccio A. V., Stinson G., Brook C. B., Wadsley J., Couch- 
man H. M. P., Shen S., Gibson B. K., Quinn T., 2012b, 
ApJL, 744, L9 

Machacek M. E., Carlson E. D., Hall L. J., 1993, in Ak- 
erlof C. W., Srednicki M. A., eds, Annals of the New 
York Academy of Sciences Vol. 688, Texas/PASCOS ’92: 
Relativistic Astrophysics and Particle Cosmology, p. 681, 
doi: 10.1111/j. 1749-6632.1993.tb43955.x 
Macias O., Gordon C., 2014, Phys. Rev. D, 89, 063515 
Madau P., Shen S., Governato F., 2014, ApJL, 789, L17 
Martizzi D., Teyssier R., Moore B., 2013, MNRAS, 432, 
1947 

Martizzi D., Faucher-Giguere C.-A., Quataert E., 2015, 
MNRAS, 450, 504 

Mashchenko S., Couchman H. M. P., Wadsley J., 2006, 
Nature, 442, 539 

Miralda-Escude J., 2002, ApJ, 564, 60 
Moore B., 1994, Nature, 370, 629 

Muratov A. L., Keres D., Faucher-Giguere C.-A., Hop¬ 
kins P. F., Quataert E., Murray N., 2015, preprint, 
(arXiv:1501.03155) 

Navarro J. F., Eke V. R., Frenk C. S., 1996, MNRAS, 283, 
L72 

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 
493 

Nelson D., Vogelsberger M., Genel S., Sijacki D., Keres D., 
Springel V., Hernquist L., 2013, MNRAS, 429, 3353 
Newman A., Treu T., Ellis R.. S., Sand D. J., Nipoti C., 
Richard J., Julio E., 2013a, in Probes of Dark Matter on 
Galaxy Scales, p. 20302 

Newman A. B., Treu T., Ellis R. S., Sand D. J., Nipoti C., 
Richard J., Julio E., 2013b, ApJ, 765, 24 
Oiiorbe J., Garrison-Kimmel S., Mailer A. H., Bullock J. S., 
Rocha M., Hahn O., 2014, MNRAS, 437, 1894 
Oiiorbe J., Boylan-Kolchin M., Bullock J. S., Hopkins P. F., 
Keres D., Faucher-Giguere C.-A., Quataert E., Murray N., 
2015, preprint, (arXiv: 1502.02036) 

Ogiya G., Mori M., 2011, ApJL, 736, L2 
Ogiya G., Mori M., 2014, ApJ, 793, 46 
Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kennicutt 
Jr. R. C., 2011, AJ, 141, 193 
Oh S.-H., et al., 2015, AJ, 149, 180 

Oppenheimer B. D., Dave R., Keres D., Fardal M., Katz 

N. , Kollmeier J. A., Weinberg D. H., 2010, MNRAS, 406, 
2325 

Papastergis E., Giovanelli R., Haynes M. P., Shankar F., 
2015, A& A, 574, A113 


© 0000 RAS, MNRAS 000, 000-000 


22 T. K. Chan et al. 


Penarrubia J., Pontzen A., Walker M. G., Koposov S. E., 

2012, ApJL, 759, L42 

Peirani S., Kay S., Silk J., 2008, A& A, 479, 123 
Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 

2013, MNRAS 430, 105 

Pontzen A., Governato F., 2012, MNRAS, 421, 3464 
Pontzen A., Governato F., 2014, Nature, 506, 171 
Porter D. H., 1985, PhD thesis, California Univ., Berkeley. 
Power C., Navarro J. F., Jenkins A., Frenk C. S., White 
S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 
338, 14 

Price D. J., 2008, Journal of Computational Physics, 227, 
10040 

Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347 
Rahmati A., Pawlik A. H., Raicevic M., Schaye J., 2013, 
MNRAS, 430, 2427 

Read J. I., Gilmore G., 2005, MNRAS, 356, 107 
Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., 
Garrison-Kimmel S., Onorbe J., Moustakas L. A., 2013, 
MNRAS, 430, 81 

Romano-Diaz E., Shlosman I., Hoffman Y., Heller C., 2008, 
ApJL, 685, L105 

Saitoh T. R., Makino J., 2013, ApJ, 768, 44 
Salucci P., Burkert A., 2000, ApJL, 537, L9 
Sijacki D., Vogelsberger M., Keres D., Springel V., Hern- 
quist L., 2012, MNRAS, 424, 2999 
Sommer-Larsen J., 2006, ApJL, 644, LI 
Spekkens K., Giovanelli R., Haynes M. P., 2005, AJ, 129, 
2119 

Spergel D. N., Bean R., Dore O., et al., 2007, ApJS, 170, 
377 

Springel V., 2005, MNRAS, 364, 1105 
Springel V., et al., 2005, Nature, 435, 629 
Steinmetz M., Navarro J., 2000, in Combes F., Mamon 
G. A., Charmandaris V., eds, Astronomical Society of the 
Pacific Conference Series Vol. 197, Dynamics of Galax¬ 
ies: from the Early Universe to the Present, p. 165 
(arXiv:astro-ph/9910573) 

Strigari L. E., Kaplinghat M., Bullock J. S., 2007, Phys. 
Rev. D, 75, 061303 

Swaters R. A., Madore B. F., van den Bosch F. C., Balcells 
M., 2003, ApJ, 583, 732 

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MN¬ 
RAS, 429, 3068 

Toilet E., et al., 2015, preprint, (arXiv: 1507.03590) 
Tonini C., Lapi A., Salucci P., 2006, ApJ, 649, 591 
Tully R. B., Fisher J. R., 1976, in Bulletin of the American 
Astronomical Society, p. 555 

Velliscig M., van Daalen M. P., Schaye J., McCarthy I. G., 
Cacciato M., Le Brun A. M. C., Dalla Vecchia C., 2014, 
MNRAS, 442, 2641 

Verheijen M. A. W., 1997, PhD thesis, PhD thesis, 
Univ. Groningen, The Netherlands , (1997) 

Vogelsberger M., Zavala J., Loeb A., 2012a, MNRAS, 423, 
3740 

Vogelsberger M., Sijacki D., Keres D., Springel V., Hern- 
quist L., 2012b, MNRAS, 425, 3024 
Walker M. G., Mateo M., Olszewski E. W., Penarrubia J., 
Wyn Evans N., Gilmore G., 2009, ApJ, 704, 1274 
Walter F., Brinks E., de Blok W. J. G., Bigiel F., Kennicutt 
Jr. R. C., Thornley M. D., Leroy A., 2008, AJ, 136, 2563 
Weisz D. R., Dolphin A. E., Skillman E. D., Holtzman J., 



Figure Al. The virial mass of a halo as a function of the max¬ 
imum particle mass allowed in order to reach convergent profiles 
at 0.3, 1, 2, 3%J? v i r , 300 pc and 700 pc. Small circles represent the 
corresponding values from dark matter only simulation listed in 
Table [I] (Our full physics simulations have the same DM particle 
numbers.) 


Gilbert K. M., Dalcanton J. J., Williams B. F., 2014, ApJ, 
789, 147 

Wheeler C., Onorbe J., Bullock J. S., Boylan-Kolchin M., 
Elbert O. D., Garrison-Kimmel S., Hopkins P. F., Keres 
D., 2015, preprint, (arXiv: 1504.02466) 

Wolf J., Martinez G. D., Bullock J. S., Kaplinghat M., 
Geha M., Munoz R.. R.., Simon J. D., Avedo F. F., 2010, 
MNRAS, 406, 1220 

Yoshida N., Springel V., White S. D. M., Tormen G., 2000, 
ApJL, 535, L103 

Zolotov A., et al., 2012, ApJ, 761, 71 
de Blok W. J. G., McGaugh S. S., 1997, MNRAS, 290, 533 
de Blok W. J. G., Walter F., Brinks E., Trachternach C., 
Oh S.-H., Kennicutt Jr. R. C., 2008, AJ, 136, 2648 
de Laix A. A., Scherrer R. J., Schaefer R. K., 1995, ApJ, 
452, 495 


APPENDIX A: CONVERGENCE LIMITS 


We use the Power convergence criterion ( Power et al.|2 003) 
to derive empirical formulae for the minimum particle mass 
needed to quantify cusps down to 0.2 — 3%R v i r . We con¬ 
sider DM only simulations and assume their profiles can be 
fitted with an NFW profile. The ratio between the scale ra¬ 
dius r s and the virial radius R v i r is determined through a 
concentration-mass relation from Dutton fe Maccid| ( 20141, 
and we use the virial overdensity from |Bryan fe Norman| 
(19981. Then we calculate the enclosed number of particles 
and density within 0.3 — 3%f? v ir as well as 300 and 700 pc. 
From Eq. [I] we estimate the minimum radius such that 
frelax/to < 0.6 and plot the relation between the required 
particle mass to meet this criteria at the desired radius, and 
the total halo mass in Figure [ATI 

The data can be fitted with a linear relation, 


logi O ("idm/M 0 ) = alog 10 (M h /M o ) + b 
where a and b are list in Table [AH 


(Al) 
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Tres / -^vir 

a 

b 

0.3% 

0.94 

-6.2 

1% 

0.95 

-5.0 

2% 

0.95 

-4.3 

3% 

0.96 

-3.9 


Table Al. The coefficients in Eq. for different resolutions. 


Oh et al. (2011) measured a between 300-700 pc in 
dwarf galaxies with M* between 10 6 — 10 9 Mq and found 
cored profiles in those galaxies. In order to match this ob¬ 
servational result, the minimum resolved radius in the simu¬ 
lations should be around 300 pc. For a fixed physical radius 
this turns into a requirement in particle mass that is almost 
independent of the halo mass, owing to higher concentration 
in lower mass halos.The maximum mass of dark matter par¬ 
ticles needs to be slightly smaller than 10 4 Mq to converge at 
300 pc and smaller than 10 5 Mq to converge within 700 pc. 
Our m09 and mlO are clearly converged at both of these 
radii. 

Our slightly more massive halos ml0hl297, 
ml0hll46 and ml0h573 are marginally converged 
at 300 pc, but fully converged at 700 pc. Figure [Xl] assumes 
an NFW profile and implicitly assumes that the central 
region has close to a Hubble time to undergo relaxation 
processes. As discussed in the main text, it is not clear 
what the appropriate convergence criterion should be once 
large cores are formed and the central density is reduced. 
This likely depends on the core formation time as well as 
details of the gravitational softening of multiple particle 
species and their time-stepping algorithms. 

It is important to note that our DM force softening 
is typically a factor of five smaller than the correspond¬ 
ing Power convergence criterion. Furthermore, we have also 
tested if the force softening of the baryonic component in¬ 
fluences dark matter profiles: e.g. we increased the baryonic 
softening from 2.0 to 25 pc in the slightly modified version 
of mlO run in Onorbe et al. (20151 and found that the dark 


matter profile was only mildly changed (the core size was 
actually larger in run with smaller softening). While two- 
body relaxation effects are important in estimating central 
DM profiles, in H14 (Appendix C) we have used idealized 
runs to show that our standard resolution in ml2 runs is 
also sufficient to reliably determine other relevant quantities, 
e.g. SFR, wind mass-loading and gas phase distribution. All 
are consistent to within a factor of ~ 2 even with ~ 50 
times better mass resolution. This indicates that the gen¬ 
eral properties of dark matter profiles on resolved scales in 
our simulations are numerically robust. 




lo glo( M h/ M Q ) 


Figure Bl. DM profile slope a inferred from different criteria 
plotted as a function of the halo mass. “0.3-0.7kpc” is the dark 
matter density slope interpolated within 0.3-0.7 kpc from the cen¬ 
ter. “1 — 2%i? v ; r ” is a interpolated within 1 — 2%R viv . “l-2kpc” 
is a interpolated within 1-2 kpc. Filled circles show that the pro¬ 
file measurement range is larger than 0.5rp ow and smaller than a 
third of r s . Open circles indicate that one of these criteria is not 
satisfied. 


meaningful, so we use this fitting range as our “default” 
choice in the main text. Overall all of the methods show very 
similar trends of the DM density profile slope with mass. 

This paper has been typeset from a TgK/ DI^HjX file prepared 
by the author. 


APPENDIX B: CHOICE OF a 

We investigate the effect of different fitting ranges on a in 
this appendix. We consider three different fitting ranges, 1-2 
kpc, 0.3-0.7 kpc, and 1 — 2%R v o. Figure [Bl] shows a that 
corresponds to those ranges. In general, 0.3-0.7 kpc includes 
some overlap below the Power radius for halos with mass 
larger than 10 11 Mq but a in this range can be directly com¬ 
pared with observations of dwarf galaxies. 1-2 kpc lies out¬ 
side the central region (> 1.3r a ) in small dwarfs. 1 — 2%R V1I 
is well-resolved, lies inside the central region and physically 
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